{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "a012e739-6f31-4721-a5d1-6160f8ee1b20",
   "metadata": {},
   "source": [
    "### Cross-validation\n",
    "This notebook demonstrates an appropriate strategy for using cross-validtion to evaluate machine learning model performance.  We utilize the strategies described in [Practically Significant Method Comparison Protocols for Machine Learning in Small Molecule Drug Discovery](https://pubs.acs.org/doi/10.1021/acs.jcim.5c01609)\n",
    "\n",
    "We begin by installing tthe necessary Python libraries. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "80e49b6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "import sys\n",
    "\n",
    "if 'google.colab' in sys.modules:\n",
    "    from IPython import get_ipython\n",
    "    # only install packages if running in Google Colab\n",
    "    get_ipython().run_line_magic('pip', 'install pandas matplotlib seaborn lightgbm scikit-learn tqdm statsmodels useful_rdkit_utils')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51da63fd",
   "metadata": {},
   "source": [
    "When performing cross-validation, it's important to prevent similar molecular structures from appearing in both the training and test sets to avoid data leakage and overestimating model performance. One effective strategy is to cluster the compounds and ensure that each cluster is assigned exclusively to either the training or the test set. Below, we use utility functions from [useful_rdkit_utils](https://github.com/PatWalters/useful_rdkit_utils) to cluster the dataset using three different approaches, each grouping molecules into clusters in a different way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f3afd3b0-1c78-4940-941b-5c9f17973206",
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import useful_rdkit_utils as uru\n",
    "from sklearn.model_selection import GroupShuffleSplit, train_test_split\n",
    "from lightgbm import LGBMRegressor\n",
    "from sklearn.metrics import r2_score, mean_absolute_error\n",
    "import numpy as np\n",
    "import warnings\n",
    "from tqdm.auto import tqdm\n",
    "import seaborn as sns\n",
    "from sklearn.base import clone\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.ensemble import HistGradientBoostingRegressor\n",
    "from statsmodels.stats.multicomp import pairwise_tukeyhsd\n",
    "from scipy.stats import f_oneway"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f66a4587",
   "metadata": {},
   "source": [
    "Read the input data. In this example we will use EC50 data for [Pregnane-X Receptor (PXR)](https://pubs.acs.org/doi/10.1021/acsmedchemlett.1c00187) agoinsm. The Pregnane-X Receptor (PXR) is a nuclear receptor that serves as a critical \"xenobiotic sensor\" in the body, functioning as a master regulator of drug metabolism and clearance. Because PXR has a notoriously large and flexible binding pocket that can accommodate a wide variety of structures, it is a common \"anti-target\" in early drug discovery. Teams frequently screen compounds to identify and filter out potent PXR activators to avoid these safety liabilities later in clinical trials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "186f0790-b9e3-43b1-b0a0-f68ade9aadbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "chembl_df = pd.read_csv(\"https://raw.githubusercontent.com/PatWalters/datafiles/refs/heads/main/PXR_chembl_EC50.csv\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5d24e4e",
   "metadata": {},
   "source": [
    "Calculate RDKit descriptors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "b23f9f12-c9b5-4ea9-aac6-952e37f8c1ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "desc = uru.RDKitDescriptors()\n",
    "chembl_df['desc'] = chembl_df.OPENADMET_CANONICAL_SMILES.apply(desc.calc_smiles)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d793e91e",
   "metadata": {},
   "source": [
    "In cross-validation, it's important to prevent similar molecular structures from appearing in both the training and test sets to avoid data leakage and overestimating model performance. One effective strategy is to cluster the compounds and ensure that each cluster is assigned exclusively to either the training or the test set. Below, we use utility functions from useful_rdkit_utils to cluster the dataset using three different approaches, each grouping molecules into clusters in a different way.\n",
    "- *random* - puts each compound in its own cluster\n",
    "- *scaffold* - uses the RDKit to decompose molecules into Bemis-Murcko clusters.  Each scaffold is then a distinct cluster. \n",
    "- *butina* - uses the RDKit to cluster compounds using the Butina algorithm. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "37dc15db-480b-46c3-8630-ef8df801cd0c",
   "metadata": {},
   "outputs": [],
   "source": [
    "random_clusters = uru.get_random_clusters(chembl_df.OPENADMET_CANONICAL_SMILES)\n",
    "scaffold_clusters = uru.get_bemis_murcko_clusters(chembl_df.OPENADMET_CANONICAL_SMILES)\n",
    "butina_clusters = uru.get_butina_clusters(chembl_df.OPENADMET_CANONICAL_SMILES)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "25d2faf6",
   "metadata": {},
   "source": [
    "Build a list of clustering methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "32c9b44e-b0a0-4c90-8de0-02387222e20d",
   "metadata": {},
   "outputs": [],
   "source": [
    "cluster_list = [random_clusters, scaffold_clusters, butina_clusters]\n",
    "method_names = [\"random\",\"scaffold\",\"butina\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f405983a",
   "metadata": {},
   "source": [
    "The function below performs nested cross-validation using GroupShuffleSplit to prevent similar groups (such as molecular clusters) from being present in both training and test sets. \n",
    "By default, it runs 5 repetitions of 5-fold cross-validation (5x5). \n",
    "Inputs:\n",
    "   - A DataFrame containing the data.\n",
    "   - A scikit-learn–compatible model.\n",
    "   - The column names for descriptors/features (`X`) and activity/target (`y`).\n",
    "   - A group label array indicating cluster membership for each compound.    \n",
    "    \n",
    " This approach provides a more realistic assessment of model performance by minimizing information leakage due to structural similarity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e63ec098-ec69-4c47-adad-0a0088dea8a9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def cross_validate(\n",
    "    df,\n",
    "    model_in,\n",
    "    X_col,\n",
    "    y_col,\n",
    "    group_labels,\n",
    "    pbar_desc=\"\",\n",
    "    metric_list=None,\n",
    "    n_outer=5,\n",
    "    n_inner=5,\n",
    "    test_size=0.2,\n",
    "):\n",
    "    \"\"\"\n",
    "    Run nested cross-validation using GroupShuffleSplit to ensure that similar groups (e.g., clusters) do not leak between train/test.\n",
    "    \n",
    "    Parameters:\n",
    "        df: pd.DataFrame - The dataset.\n",
    "        model_in: estimator - The estimator to fit.\n",
    "        X_col: str - The feature column (containing numpy vectors).\n",
    "        y_col: str - The target column.\n",
    "        group_labels: list/array - Cluster or group assignment used for splitting.\n",
    "        pbar_desc: str - Optional progress bar description.\n",
    "        metric_list: list of callables - Scoring functions; defaults to [r2_score].\n",
    "        n_outer: int - Number of outer iterations (CV repeats).\n",
    "        n_inner: int - Number of splits per outer iteration.\n",
    "        test_size: float - Test fraction for GroupShuffleSplit.\n",
    "\n",
    "    Returns:\n",
    "        result_list: list - Metric results, shape [n_outer * n_inner, len(metric_list)]\n",
    "    \"\"\"\n",
    "    from sklearn.base import clone\n",
    "    import warnings\n",
    "    from tqdm.auto import tqdm\n",
    "    import numpy as np\n",
    "\n",
    "    if metric_list is None:\n",
    "        from sklearn.metrics import r2_score\n",
    "        metric_list = [r2_score]\n",
    "\n",
    "    result_list = []\n",
    "    y_arr = df[y_col].values\n",
    "    X_vals = np.stack(df[X_col])\n",
    "    for outer_idx in tqdm(range(n_outer), desc=pbar_desc):\n",
    "        gss = GroupShuffleSplit(n_splits=n_inner, test_size=test_size, random_state=None)\n",
    "        for train_idx, test_idx in gss.split(X_vals, y=y_arr, groups=group_labels):\n",
    "            # Use index on the input df\n",
    "            train = df.iloc[train_idx]\n",
    "            test = df.iloc[test_idx]\n",
    "\n",
    "            model = clone(model_in)\n",
    "            model.fit(np.stack(train[X_col]), train[y_col].values)\n",
    "\n",
    "            with warnings.catch_warnings():\n",
    "                warnings.simplefilter(\"ignore\", category=UserWarning)\n",
    "                pred = model.predict(np.stack(test[X_col]))\n",
    "\n",
    "            results_row = []\n",
    "            for metric in metric_list:\n",
    "                try:\n",
    "                    results_row.append(metric(test[y_col].values, pred))\n",
    "                except Exception as e:\n",
    "                    results_row.append(np.nan)\n",
    "                    print(f\"Warning: Metric {metric.__name__} failed: {e}\")\n",
    "            result_list.append(results_row)\n",
    "\n",
    "    return result_list"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64bd4c5c",
   "metadata": {},
   "source": [
    "Next, we will perform cross-validation using each of the three data splitting strategies. The example below uses `LGBMRegressor`, but you can substitute any regressor that is compatible with scikit-learn. We define `metric_list` to specify which evaluation metrics to calculate during validation. For each splitting strategy, the resulting scores along with their associated method names are collected and stored in the `split_score_df_list`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6040df6e-e85d-4e3f-ba5f-a59fdae5b1c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "629b2b1c2f394f969d1c6c4ad08948ba",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "random:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b45fcd1539574e59ada33291310dc7cc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "scaffold:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9b2de8d254374523a746f36617646c22",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "butina:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "current_model = LGBMRegressor(verbose=-1)\n",
    "metric_list = [r2_score,mean_absolute_error]\n",
    "metric_names = [\"R2\",\"MAE\"]\n",
    "split_score_df_list = []\n",
    "for cluster,method_name in zip(cluster_list,method_names):\n",
    "    split_score_list = cross_validate(chembl_df,current_model,\"desc\",\"PXR_pEC50\",cluster,method_name,metric_list)\n",
    "    split_score_df = pd.DataFrame(split_score_list,columns=metric_names)\n",
    "    split_score_df['method'] = method_name\n",
    "    split_score_df_list.append(split_score_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a95e2a14",
   "metadata": {},
   "source": [
    "Combine the results into a single DataFrame for easier comparison."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "7bd14090-fb5b-42da-9fcd-11669039038c",
   "metadata": {},
   "outputs": [],
   "source": [
    "split_combo_df = pd.concat(split_score_df_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "852f5af3",
   "metadata": {},
   "source": [
    "Plot $R^2$ and MAE as box plots.  Note that random splits and Butina splits tend to provide uppper and lower bounds for model performance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "91745519-1127-4892-a1aa-a7293f27a51d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAPxZJREFUeJzt3Q+cVWWdP/Avw1/xz6ASoMgf0VCRBIMgNK02DNeybFtDUyFSWP9QbpQpi0JoI2VFtIoSrmSarqZl626GurToz8SlIHetlM1/A78UhDUGRQFj5vd6zu8182JkwIGZc+/Mve/363Vewzn3nHueO/dyn/mc5znP06Gurq4uAAAAgFZX0fpPCQAAACRCNwAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA56ZTXE5eq2traeOmll2L//fePDh06FLs4AJShurq6eO211+LQQw+NioryvX6uTgagPdTHQvceSpV7v379il0MAIg1a9bEYYcdFuVKnQxAe6iPhe49lK6m1/9iDzjggGIXB4AytGnTpixs1tdJ5UqdDEB7qI+F7j1U330tVe4qeACKqdy7VKuTAWgP9XH53ggGAAAAORO6AQAAICdCNwAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA5EboBAAAgJ0I3AAAA5EToBgAAgJwI3QBAq5g/f34MHDgwunXrFqNHj47ly5fvdv958+bFUUcdFfvss0/069cvvvSlL8WWLVsKVl4AKAShGwBosbvvvjumTZsWs2bNipUrV8awYcNi3Lhx8corrzS5/5133hlXXHFFtv/TTz8dt9xyS/Yc//AP/1DwsgNAnoRuAKDF5s6dG5MnT45JkybFkCFDYsGCBdG9e/dYtGhRk/s//vjjceKJJ8ZnP/vZrHX8ox/9aJx99tnv2DoOAO1Np2IXgHyk7nmrV6+OUte/f/+sGyMAxbNt27ZYsWJFTJ8+vWFbRUVFjB07NpYtW9bkMSeccEL86Ec/ykL2qFGj4vnnn48HHnggzjvvvF2eZ+vWrdlSb9OmTa38SgDaB3/rty9Cd4lK/wmnTJkSpW7hwoUxePDgYhcDoKxt2LAhtm/fHr179260Pa0/88wzTR6TWrjTcR/4wAeirq4u/vKXv8SFF1642+7lc+bMidmzZ7d6+QHaG3/rty9CdwlfFUof0kKprq6OqqqqmDFjRgwYMKCgrxOA9mfp0qVx7bXXxo033pgNuvbss8/GpZdeGtdcc01cddVVTR6TWtLTfeM7tnSnAdgAyo2/9dsXobtEpW4YxbgqlP4TlsLVKACar2fPntGxY8dYt25do+1pvU+fPk0ek4J16kp+wQUXZOvvec97YvPmzVnLTfqjLnVPf7uuXbtmC0C587d++2IgNQCgRbp06RIjRoyIJUuWNGyrra3N1seMGdPkMW+88cZOwToF9yR1NweAUqGlGwBosdTte+LEiTFy5MhsYLQ0B3dquU6jmScTJkyIvn37ZvdlJ6effno24vnxxx/f0L08tX6n7fXhGwBKgdANALTY+PHjY/369TFz5sxYu3ZtDB8+PBYvXtwwuFoa9GfHlu0rr7wyOnTokP3805/+FO9617uywJ3uGQSAUiJ0AwCtYurUqdmyq4HTdtSpU6eYNWtWtgBAKXNPNwAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA5EboBAAAgJ0I3AAAA5EToBgAAgJwI3QAAAJAToRsAAAByInQDAABAToRuAAAAyInQDQAAADkRugEAACAnQjcAAADkpCRC9/z582PgwIHRrVu3GD16dCxfvny3+2/cuDEuueSSOOSQQ6Jr164xePDgeOCBBwpWXgAAAMpDp2jn7r777pg2bVosWLAgC9zz5s2LcePGxapVq6JXr1477b9t27Y45ZRTssfuvffe6Nu3b1RXV0ePHj2KUn4AAABKV7sP3XPnzo3JkyfHpEmTsvUUvn/+85/HokWL4oorrthp/7T91Vdfjccffzw6d+6cbUut5AAAANDa2nX38tRqvWLFihg7dmzDtoqKimx92bJlTR5z//33x5gxY7Lu5b17946hQ4fGtddeG9u3b29y/61bt8amTZsaLQAAAFDyLd0bNmzIwnIKzztK688880yTxzz//PPxy1/+Ms4555zsPu5nn302Lr744njrrbdi1qxZO+0/Z86cmD17dquUd926dVFTUxOlKHXR3/FnKaqsrNzpswYAAFCyoXtv1NbWZvdzL1y4MDp27BgjRoyIP/3pT/Gtb32rydA9ffr07J7xeqmlu1+/fnsVuM89b0K8tW1rlLKqqqooVZ27dI0f3X6b4A0AAJRH6O7Zs2cWnFOg3VFa79OnT5PHpBHL073c6bh6xxxzTKxduzbrrt6lS5dG+6fRzdPSUqmFOwXuNwd9MGq7Vbb4+Sisii01Ec8/kr2PQjcAAFAWoTsF5NRSvWTJkjjjjDMaWrLT+tSpU5s85sQTT4w777wz2y/d/538z//8TxbG3x6485ACd+2+PXM/DwAAAMXXrgdSS1LX75tvvjl++MMfxtNPPx0XXXRRbN68uWE08wkTJmRdxOulx9Po5ZdeemkWttNI52kgtTSwGgAAALSmdt3SnYwfPz7Wr18fM2fOzLqIDx8+PBYvXtzQBXj16tUNLdpJuh/7wQcfjC996Utx3HHHZfN0pwB++eWXF/FVAAAAUIrafehOUlfyXXUnX7p06U7b0pRhTzzxRAFKBgAAQDlr993LAQAAoK0qiZZuoGW2bNmS3YpRyvr37x/dunUrdjEAACgzQjeQBe4pU6ZEKVu4cGEMHjy42MUAAKDMCN1A1gqcQmkhVFdXR1VVVcyYMSMGDBgQhXyNAABQaEI3kHW7LnQrcArcWp4BACh1BlIDAACAnAjdAAAAkBOhGwAAAHLinu4Cq3hzY7GLwF7wvgEAAHtD6C6wfV54tNhFAAAAoECE7gJ78/CTo3afHsUuBnvR0u2CCQAAsKeE7gJLgbt2357FLgYAAAAFYCA1AAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnAjdAAAAkBOhGwAAAHIidAMAAEBOhG4AAADIidANAAAAORG6AQAAICdCNwAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA5EboBAAAgJ53yemKgZdatWxc1NTVRaqqrqxv9LEWVlZXRu3fvYhcDAIA2QOiGNhq4zz1vQry1bWuUqqqqqihVnbt0jR/dfpvgDQCA0A1tUWrhToH7zUEfjNpulcUuDnugYktNxPOPZO+h0A0AgNBdjD/IaXeK9b6lwF27b8+inBsAAGg5obuA93imLqepBYz2Kb1/6X0EAABoLqG7QFI303SPZykOjFU/KFa6R3fGjBkxYMCAKEUGxwIAAPaU0F1AKbCVemhLgXvw4MHFLgYAAECbYJ5uAAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnAjdAAAAkBOhGwAAAHIidAMAAEBOhG4AoFXMnz8/Bg4cGN26dYvRo0fH8uXLd7nvhz70oejQocNOy8c+9rGClhkA8iZ0AwAtdvfdd8e0adNi1qxZsXLlyhg2bFiMGzcuXnnllSb3/+lPfxovv/xyw/K73/0uOnbsGGeeeWbByw4AeRK6AYAWmzt3bkyePDkmTZoUQ4YMiQULFkT37t1j0aJFTe5/0EEHRZ8+fRqWhx9+ONtf6Aag1AjdAECLbNu2LVasWBFjx45t2FZRUZGtL1u2rFnPccstt8RZZ50V++677y732bp1a2zatKnRAgBtndANALTIhg0bYvv27dG7d+9G29P62rVr3/H4dO936l5+wQUX7Ha/OXPmRGVlZcPSr1+/FpcdAPImdAMARZVaud/znvfEqFGjdrvf9OnTo6ampmFZs2ZNwcoIAHur014fCQAQET179swGQVu3bl2j7Wk93a+9O5s3b4677rorrr766nc8T9euXbMFANoTLd0AQIt06dIlRowYEUuWLGnYVltbm62PGTNmt8fec8892b3a5557bgFKCgCFp6UbAGixNF3YxIkTY+TIkVk38Xnz5mWt2Gk082TChAnRt2/f7L7st3ctP+OMM+Lggw8uUskBIF9CNwDQYuPHj4/169fHzJkzs8HThg8fHosXL24YXG316tXZiOY7WrVqVTz22GPx0EMPFanUAJA/oRsAaBVTp07NlqYsXbp0p21HHXVU1NXVFaBkAFA87ukGAACAnAjdAAAAkBOhGwAAAHIidAMAAEBOhG4AAADIidHLAWi2LVu2ZFM/lbr+/ftHt27dil0MAKAECN0ANFsK3FOmTIlSt3Dhwhg8eHCxiwEAlAChG4A9agFOgbRQqquro6qqKmbMmBEDBgwo6OsEAGgNQjcAzZa6XBejBTgFbi3PAEB7ZCA1AAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnAjdAAAAkJOSCN3z58+PgQMHZlPZjB49OpYvX96s4+66667o0KFDnHHGGbmXEQAAgPLT7kP33XffHdOmTYtZs2bFypUrY9iwYTFu3Lh45ZVXdnvciy++GF/5ylfipJNOKlhZAQAAKC/tPnTPnTs3Jk+eHJMmTYohQ4bEggULonv37rFo0aJdHrN9+/Y455xzYvbs2TFo0KCClhcAAIDy0a5D97Zt22LFihUxduzYhm0VFRXZ+rJly3Z53NVXXx29evWK888//x3PsXXr1ti0aVOjBQAAAEo+dG/YsCFrte7du3ej7Wl97dq1TR7z2GOPxS233BI333xzs84xZ86cqKysbFj69evXKmUHAACg9LXr0L2nXnvttTjvvPOywN2zZ89mHTN9+vSoqalpWNasWZN7OQEAACgNnaIdS8G5Y8eOsW7dukbb03qfPn122v+5557LBlA7/fTTG7bV1tZmPzt16hSrVq2KI444otExXbt2zRYAAAAoq5buLl26xIgRI2LJkiWNQnRaHzNmzE77H3300fHUU0/Fk08+2bB84hOfiA9/+MPZv3UdBwAAoDW165buJE0XNnHixBg5cmSMGjUq5s2bF5s3b85GM08mTJgQffv2ze7NTvN4Dx06tNHxPXr0yH6+fTsAAABEuYfu8ePHx/r162PmzJnZ4GnDhw+PxYsXNwyutnr16mxEcwAAACi0dh+6k6lTp2ZLU5YuXbrbY2+99dacSgUAAEC50wQMAAAAOSmJlm52tmXLlqxrfaFUV1c3+lko/fv3z+7VBwAAaIuE7hKVAveUKVMKft6qqqqCnm/hwoUxePDggp4TAACguYTuEpVagFMgLYfXCQAA0FYJ3SUqdbnWAgwAAFBcBlIDAACAnAjdAAAAkBOhGwAAAHIidAMAAEBODKQGbVjFmxuLXQT2kPcMAIAdCd3Qhu3zwqPFLgIAANACQje0YW8efnLU7tOj2MVgD1u6XSwBAKCe0A1tWArctfv2LHYxAACAvWQgNQAAAMiJ0A0AAAA5EboBAAAgJ0I3AAAA5EToBgAAgJwI3QAAAJAToRsAAAByInQDAABAToRuAAAAyInQDQAAADkRugEAACAnQjcAAADkpFNeTwwAAFBO1q1bFzU1NVFqqqurG/0sRZWVldG7d+9cnlvoBgAAaIXAfe55E+KtbVujVFVVVUWp6tyla/zo9ttyCd5CN0A7V6pX1RNX1gFoL1JdnAL3m4M+GLXdKotdHPZAxZaaiOcfyd5DoRuAsruqnriyDkB7kQJ37b49i10M2hChG6Adc1W9fcv7yjoAUHxCN0AJcFUdAKBtMmUYAAAA5EToBgAAgJwI3QAAAJAToRsAAAByInQDAABAToRuAAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnAjdAAAAkBOhGwAAAHIidAMAAEBOOuX1xEDLVWypKXYR2EPeMwAAdiR0QxtUWVkZnbt0jXj+kWIXhb2Q3rv0HgIAgNANbVDv3r3jR7ffFjU1pddqWl1dHVVVVTFjxowYMGBAlKIUuNN7CAAAQje0USm0lXJwS4F78ODBxS4GAADkykBqAECrmD9/fgwcODC6desWo0ePjuXLl+92/40bN8Yll1wShxxySHTt2jW7EPfAAw8UrLwAUAhaugGAFrv77rtj2rRpsWDBgixwz5s3L8aNGxerVq2KXr167bT/tm3b4pRTTskeu/fee6Nv377Z7Sc9evQoSvkBIC9CNwDQYnPnzo3JkyfHpEmTsvUUvn/+85/HokWL4oorrthp/7T91Vdfjccffzw6d+6cbUut5ABQanQvBwBaJLVar1ixIsaOHduwraKiIltftmxZk8fcf//9MWbMmKx7eRq/YujQoXHttdfG9u3bd3merVu3xqZNmxotANDWCd0AQIts2LAhC8tvH/wxra9du7bJY55//vmsW3k6Lt3HfdVVV8V3vvOd+PrXv77L88yZMyebHaB+6devX6u/FgBobUI3AFBwtbW12f3cCxcujBEjRsT48eOzqQRTt/RdmT59ejaVYv2yZs2agpYZAPaGe7oBgBbp2bNndOzYMdatW9doe1rv06dPk8ekEcvTvdzpuHrHHHNM1jKeuqt36dJlp2PSCOdpAYD2ROimxVLXwP/+7//OBsQ56KCD4rjjjmv0RxQApS0F5NRavWTJkjjjjDMaWrLT+tSpU5s85sQTT4w777wz2y/d/538z//8TxbGmwrcUA62bNkSq1evjlLXv3//bGpBKBdCNy3y6KOPxo033tjonr3UqnHxxRfHySefXNSyAVA4abqwiRMnxsiRI2PUqFHZlGGbN29uGM18woQJ2bRg6b7s5KKLLoobbrghLr300vjCF74Qf/zjH7OB1L74xS8W+ZVA8aTAPWXKlCh16baSwYMHF7sYUDBCNy0K3LNmzcpGn00D4Bx++OHxwgsvxB133JFtnz17tuANUCbSPdnr16+PmTNnZhdihw8fHosXL24YXC2FifoW7SQNgvbggw/Gl770payHVArkKYBffvnlRXwVUPwW4BRIC6W6ujqqqqqy8RQGDBhQ0NcJ5UToZq+7lKcW7hS400iz9X9IHXvssdn6lVdeGTfddFPWfVBXc4DykLqS76o7+dKlS3faluqQJ554ogAlg/YhdbkuRgtwCtxaniE/Ri9nr6R7uFNLxjnnnNOo5SJJ62n7yy+/nO0HAABQroRu9koaNC1JXcqbUr+9fj8AAIByJHSzV9Io5Um6h7sp9dvr9wMAAChHQjd7JQ16k0YpT4OmpeledpTW0/Y07UvaDwAAoFwJ3eyVNDhamhZs2bJl2aBpv//97+ONN97Ifqb1tD1NB2MQNQAAoJwJ3ey1NB1Ymhbs+eefj0suuSROO+207GfqWm66MIC2bfny5dlMFLuydevW+PGPf1zQMgFAKTJlGC2SgnWaFiyNUp4GTUv3cKcu5Vq4Adq2NF1XmmWiV69e2foBBxwQTz75ZAwaNChb37hxY5x99tnxmc98psglBYD2TeimxVLAPv7444tdDAD2QF1d3W7Xd7UNANgzupcDAE3q0KFDsYsAAO2e0A0AAABtpXv5m2++md2727dv30bb06jVxx57bGuWDQDI0R/+8IdYu3ZtQ1fyZ555Jl5//fVsfcOGDUUuHQCUYei+99574+///u+jZ8+e2VzMN998c4wePTp77LzzzouVK1fmVU4AoJV95CMfaXTf9sc//vGGbuVpu+7lAFDg7uVf//rXY8WKFdnopj/4wQ/i/PPPjzvvvLPog63Mnz8/Bg4cGN26dcsuAqRpUHYlXSg46aST4sADD8yWsWPH7nZ/AChFaXrHNOVj+vn2pX57+gkAFLCl+6233orevXtn/x4xYkQ8+uij8alPfSqeffbZol0Nv/vuu2PatGmxYMGCLHDPmzcvxo0bF6tWrWqYBmVHS5cuzaZAOeGEE7KQ/s1vfjM++tGPZt3j395lHgBK1YABA95xn9/97ncFKQsAlLI9aulOITbNx1wvzcn88MMPx9NPP91oeyHNnTs3Jk+eHJMmTYohQ4Zk4bt79+6xaNGiJve/44474uKLL47hw4fH0UcfHf/0T/+UdZVfsmRJwcsOAG3Na6+9FgsXLoxRo0bFsGHDil0cACiv0H377bfv1HrcpUuX+Od//ud45JFHotC2bduWdXdPXcTrVVRUZOvLli1r1nO88cYbWQt+uoAAAOUq9V6bOHFiHHLIIfHtb387/uqv/iqeeOKJYhcLAMqre/lhhx22y8dOPPHEKLQ0sur27dsburzXS+tpBNbmuPzyy+PQQw9tFNx3tHXr1mypt2nTphaWGgDahjRy+a233hq33HJLVr995jOfyeq8n/3sZ1nvMQCgyPN0V1dXx0MPPdQw3cjbvfTSS9GWfeMb34i77ror7rvvvuz+7qbMmTMnKisrG5Z+/foVvJwA0NpOP/30OOqoo7Lbw9J4KKnOvv7664tdLAAoOXsdulOX8iOPPDJOPfXUGDRoUNb1PFm9enUWZtOgZv379488panLOnbsGOvWrWu0Pa336dNnt8emrnOpnOmiwXHHHbfL/aZPnx41NTUNy5o1a1qt/ABQLL/4xS+yWUhmz54dH/vYx7L6FABoQ6H7mmuuiS984Qvx1FNPxSmnnBIXXXRRXHXVVXHEEUdkXdVGjhwZ99xzT+Qp3U+eRlHfcRC0+kHRxowZs8vjrrvuuqz8ixcvzsq5O127do0DDjig0QIA7d1jjz2WDZqW6tF0ofyGG27IbtsCAIp4T/eOnnvuubj00kuzKUfSPNmpVftXv/pV1k3tmGOOiUJJ04WlgV9SeE4jraYucps3b85GM08mTJiQTQWWuoknaYqwmTNnZvOLp7m967vG77ffftkCAOXg/e9/f7akejNNv5lm/Uh1arp4nWYmSbdT7b///sUuJgCUb0t3GvF7n332aRhgLd0TnbpsFzJwJ+PHj8/Om4J0mgbsySefzFqw6wdXS93dX3755Yb9b7rppmzU87/927/NRmitX9JzAEC52XfffePzn/981vKdeq99+ctfzm6/SrOVfOITnyh28QCgfFu6k9RanO7pTvNdp3vBDjzwwCiGqVOnZktTli5d2mj9xRdfLFCpAKB9SQOrpVuwUu+wf/u3f8tavwGAIoXuk046KWbNmpVdEU9he8uWLfG9730vTjjhhBg6dGgMHjw4OnVqUaYHAHKSWrffycEHH1yQsgBAKdvrVPzII49kP//4xz/GihUrYuXKldly2223xcaNG7NBzlLwTvd4AwBtSxr0NI3Lcvzxx0ddXV2T+3To0KHg5QKAUtPipuh3v/vd2XLWWWc1bHvhhRfiN7/5Tfz2t79t6dMDADlIs46k6T9TnZ0GHz333HPjoIMOKnaxAKDk7PVAartz+OGHx5lnnhnXXnttHk8PALRQmnkkDTT61a9+Nf71X/81G638M5/5TDz44IO7bPkGANpI6AYA2r6uXbvG2WefnU0R9oc//CGOPfbYuPjii7MpNV9//fViFw8ASoLQDQBERUVFdg93auXevn17sYsDACVD6AaAMrV169bsvu5TTjklG/w0zdN9ww03xOrVq2O//fYrdvEAoCSY0wsAylDqRn7XXXdl93Kn6cNS+O7Zs2exiwUAJUfoBoAytGDBgujfv38MGjQomwa0firQt/vpT39a8LIBQCkRugGgDE2YMME83ABQAEI3AJShW2+9tdhFAICyYCA1AAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnAjdAAAAkBOhGwAAAHJinm4AAIBWUvHmxmIXgTb2ngndQGzZsiVWr15dkHNVV1c3+lko/fv3j27duhX0nABA+dnnhUeLXQTaGKEbyAL3lClTCnrOqqqqgp5v4cKFMXjw4IKeEwAoP28efnLU7tOj2MVgD1u687xYInQDWStwCqWl/hoBAPKWAnftvj2LXQzaEKEbyLpdawUGAIDWZ/RyAAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnAjdAAAAkBOhGwAAAHIidAMAAEBOhG4AAADISae8nhgAANq7devWRU1NTZSi6urqRj9LUWVlZfTu3bvYxaDMCd0AALCLwH3ueRPirW1bo5RVVVVFqercpWv86PbbBG+KSugGAIAmpBbuFLjfHPTBqO1WWezisIcqttREPP9I9j4K3RST0A0AALuRAnftvj2LXQygnTKQGgAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA5EboBAAAgJ0I3AAAA5EToBgAAgJx0yuuJAQAAyk3FlppiF4E29p4J3QAAAC1UWVkZnbt0jXj+kWIXhb2Q3rv0HuZB6AYAAGih3r17x49uvy1qakqvpbu6ujqqqqpixowZMWDAgChFlZWV2XuYB6EbACBHW7ZsidWrV0cp69+/f3Tr1q3YxYCiS6Etr+DWFqTAPXjw4GIXo90RugEAcpQC95QpU6KULVy40B/iALsgdAMArWL+/PnxrW99K9auXRvDhg2L66+/PkaNGtXkvrfeemtMmjSp0bauXbtmrcKl2AqcQmkpdwNNrxGApgndAECL3X333TFt2rRYsGBBjB49OubNmxfjxo2LVatWRa9evZo85oADDsger9ehQ4coRanbdTFagXUDBWgbzNMNALTY3LlzY/LkyVnr9ZAhQ7Lw3b1791i0aNEuj0khu0+fPg1LKd8HCUD50tINUAIq3txY7CJQxu/btm3bYsWKFTF9+vSGbRUVFTF27NhYtmzZLo97/fXXs9bY2traeO973xvXXnttHHvssbvcf+vWrdlSb9OmTa34KgAgH0I3QAnY54VHi10EytiGDRti+/btO7VUp/VnnnmmyWOOOuqorBX8uOOOy6bX+fa3vx0nnHBC/P73v4/DDjusyWPmzJkTs2fPzuU1AEBehG6AEvDm4SdH7T49il0M9qKlu1wvmIwZMyZb6qXAfcwxx8T3v//9uOaaa5o8JrWkp/vGd2zp7tevX0HKCwB7S+gGKAEpcNfu27PYxaBM9ezZMzp27Bjr1q1rtD2tp3u1m6Nz585x/PHHx7PPPrvLfdLo5mkBgPbEQGoAQIt06dIlRowYEUuWLGnYlu7TTus7tmbvTuqe/tRTT8UhhxySY0kBoPC0dAMALZa6fU+cODFGjhyZzc2dpgzbvHlzw1zcEyZMiL59+2b3ZSdXX311vP/9748jjzwyNm7cmM3vneaXvuCCC4r8SgCgdQndAECLjR8/PtavXx8zZ86MtWvXxvDhw2Px4sUNg6utXr06G9G83p///OdsirG074EHHpi1lD/++OPZdGMAUEqEbgCgVUydOjVbmrJ06dJG69/97nezBQBKnXu6AQAAICdCNwAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA5EboBAAAgJ53yemIAACgFFW9uLHYR2AveN9qKkgjd8+fPj29961uxdu3aGDZsWFx//fUxatSoXe5/zz33xFVXXRUvvvhivPvd745vfvObcdpppxW0zAAAtA/7vPBosYsAtGPtPnTffffdMW3atFiwYEGMHj065s2bF+PGjYtVq1ZFr169dtr/8ccfj7PPPjvmzJkTH//4x+POO++MM844I1auXBlDhw4tymsAAKDtevPwk6N2nx7FLgZ70dLtggltQbsP3XPnzo3JkyfHpEmTsvUUvn/+85/HokWL4oorrthp/+9973tx6qmnxmWXXZatX3PNNfHwww/HDTfckB0LAAA7SoG7dt+exS4G0E6164HUtm3bFitWrIixY8c2bKuoqMjWly1b1uQxafuO+yepZXxX+wMAAEBZtnRv2LAhtm/fHr179260Pa0/88wzTR6T7vtuav+0vSlbt27NlnqbNm1qlbIDAABQ+tp1S3chpHu/KysrG5Z+/foVu0gAAAC0E+06dPfs2TM6duwY69ata7Q9rffp06fJY9L2Pdl/+vTpUVNT07CsWbOmFV8BAAAApaxdh+4uXbrEiBEjYsmSJQ3bamtrs/UxY8Y0eUzavuP+SRpIbVf7d+3aNQ444IBGCwAAAJT8Pd1Jmi5s4sSJMXLkyGxu7jRl2ObNmxtGM58wYUL07ds36yaeXHrppfHBD34wvvOd78THPvaxuOuuu+I3v/lNLFy4sMivBAAAgFLT7kP3+PHjY/369TFz5sxsMLThw4fH4sWLGwZLW716dTaieb0TTjghm5v7yiuvjH/4h3+Id7/73fGzn/3MHN0AAAC0unYfupOpU6dmS1OWLl2607YzzzwzWwAAACBP7fqebgAAAGjLhG4AAADIidANAAAAORG6AQAAICdCNwAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA56ZTXEwMAtFXr1q2LmpqaKEXV1dWNfpaaysrK6N27d0HPWbGlND8rpc77RlshdAMAZRe4zz1vQry1bWuUsqqqqihFnbt0jR/dfltBgncK+Ol88fwjuZ+LfKT3L72PUExCNwBQVlILdwrcbw76YNR288d4u2u5fP6R7D0sROhO50gBv5R7RaSLMzNmzIgBAwZEKSpGzwh4O6EbAChLKXDX7tuz2MWgjUuBrdRDWwrcgwcPLnYxoGQZSA0AAAByInQDAABAToRuAAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnAjdAAAAkJNOeT0xAEBbVvHmxmIXgT3kPQPaI6EbAChL+7zwaLGLAEAZELoBgLL05uEnR+0+PYpdDPawpdvFEqC9EboBgLKUAnftvj2LXQwASpyB1AAAACAnQjcAAADkROgGAACAnAjdAAAAkBOhGwAAAHIidAMAAEBOhG4AAADIidANAAAAORG6AQAAICed8npiAAqnYktNsYvAXvC+AUDpE7oB2rHKysro3KVrxPOPFLso7KX0/qX3EQAoTUI3QDvWu3fv+NHtt0VNTWm2mFZXV0dVVVXMmDEjBgwYEKUoBe70PgIApUnoBmjnUmAr9dCWAvfgwYOLXQyAXG3ZsiVWr15d0AubO/4slP79+0e3bt0Kek4oJqEbAADagBS4p0yZUvDzph5FhbRw4UIXUikrQjcAALQBqQU4BdJyeJ1QToRuAABoA1KXay3AUHrM0w0AAAA5EboBAAAgJ0I3AAAA5EToBgAAgJwYSA0AKEsVW2qKXQT2kPcMaI+EbgCgrFRWVkbnLl0jnn+k2EVhL6T3Lr2HAO2F0A0AlJXevXvHj26/LWpqSrPVtLq6OqqqqmLGjBkxYMCAKDUpcKf3EKC9ELoBgLKTQlupB7cUuM35DFB8BlIDAACAnAjdAAAAkBOhGwBoFfPnz4+BAwdGt27dYvTo0bF8+fJmHXfXXXdFhw4d4owzzsi9jABQaEI3ANBid999d0ybNi1mzZoVK1eujGHDhsW4cePilVde2e1xL774YnzlK1+Jk046qWBlBYBCEroBgBabO3duTJ48OSZNmhRDhgyJBQsWRPfu3WPRokW7PGb79u1xzjnnxOzZs2PQoEEFLS8AFIrQDQC0yLZt22LFihUxduzYhm0VFRXZ+rJly3Z53NVXXx29evWK888/v1nn2bp1a2zatKnRAgBtndANALTIhg0bslbrt0/BldbXrl3b5DGPPfZY3HLLLXHzzTc3+zxz5szJ5miuX/r169fisgNA3oRuAKCgXnvttTjvvPOywN2zZ89mHzd9+vSoqalpWNasWZNrOQGgNXRqlWcBAMpWCs4dO3aMdevWNdqe1vv06bPT/s8991w2gNrpp5/esK22tjb72alTp1i1alUcccQROx3XtWvXbAGA9kRLNwDQIl26dIkRI0bEkiVLGoXotD5mzJid9j/66KPjqaeeiieffLJh+cQnPhEf/vCHs3/rNg5AKdHSDQC0WJoubOLEiTFy5MgYNWpUzJs3LzZv3pyNZp5MmDAh+vbtm92XnebxHjp0aKPje/Tokf18+3YAaO+EbgCgxcaPHx/r16+PmTNnZoOnDR8+PBYvXtwwuNrq1auzEc0BoNwI3QBAq5g6dWq2NGXp0qW7PfbWW2+NUrVly5bsokOhVFdXN/pZCP379896MACwM6EbACBHKXBPmTKl4Oetqqoq2LkWLlwYgwcPLtj5ANoToRsAIOdW4BRKS/01AtA0oRsAIEep27VWYIDyZUQTAAAAyInQDQAAADkRugEAACAnQjcAAADkROgGAACAnLTr0ctfffXV+MIXvhD/+q//GhUVFfHpT386vve978V+++23y/1nzZoVDz30UDZn5rve9a4444wz4pprronKysqClx8AAGBPbdmyJcszhVJdXd3oZyGnI+zWrVu0d+06dJ9zzjnx8ssvx8MPPxxvvfVWTJo0KaZMmRJ33nlnk/u/9NJL2fLtb387hgwZkn1oLrzwwmzbvffeW/DyAwAA7KkUuFPuKbSqqqqCnm/hwoUlMeViuw3dTz/9dCxevDh+/etfx8iRI7Nt119/fZx22mlZqD700EN3Ombo0KHxk5/8pGH9iCOOyD445557bvzlL3+JTp3a7a8DAAAoE6kFOAXScnidpaDdpsxly5ZFjx49GgJ3Mnbs2Kyb+X/+53/Gpz71qWY9T01NTRxwwAECNwAA0C6kLtel0AJcLtpt0ly7dm306tWr0bYUnA866KDssebYsGFDdj/37rpmbN26NVvqbdq0qQWlBgAAoJy0udHLr7jiiujQocNul2eeeabF50nh+WMf+1h2b/fXvva1Xe43Z86cbJC1+qVfv34tPjcAAADloc21dH/5y1+Oz33uc7vdZ9CgQdGnT5945ZVXGm1P92WnEcrTY7vz2muvxamnnhr7779/3HfffdG5c+dd7jt9+vSYNm1ao7AueAMAANAuQ3eaxist72TMmDGxcePGWLFiRYwYMSLb9stf/jJqa2tj9OjRuzwuheZx48ZF165d4/7773/HIejTfmkBAACAdt+9vLmOOeaYrLV68uTJsXz58vjVr34VU6dOjbPOOqth5PI//elPcfTRR2eP1wfuj370o7F58+a45ZZbsvV0/3datm/fXuRXBAAAQKlpcy3de+KOO+7IgvZHPvKRbNTyT3/60/GP//iPDY+nubtXrVoVb7zxRra+cuXKbGTz5Mgjj2z0XC+88EIMHDiwwK8AAACAUtauQ3caqfzOO+/c5eMpRNfV1TWsf+hDH2q0DgAAAHlqt93LAQAAoK0TugEAACAnQjcAAADkROgGAACAnAjdAAAAkBOhGwAAAHIidAMAAEBOhG4AAADIidANAAAAORG6AQAAICdCNwAAAORE6AYAAICcCN0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA56ZTXEwMAUFjbt2+P//7v/45XX301DjrooDjuuOOiY8eOxS4WQFkTugEASsCjjz4aN954Y6xdu7ZhW58+feLiiy+Ok08+uahlAyhnupcDAJRA4J41a1YMGjQo5s+fHw888ED2M62n7elxAIpD6AYAaOddylML95gxY+LrX/96HHvssdG9e/fsZ1pP22+66aZsPwAKT/dyAJpty5YtsXr16oKdr7q6utHPQunfv39069atoOeEvZXu4U5dyq+66qqoqGjcnpLWzznnnLjkkkuy/Y4//viilROgXAndADRbCtxTpkwp+HmrqqoKer6FCxfG4MGDC3pO2Ftp0LTk8MMPb/Lx+u31+wFQWEI3AHvUApwCaTm8Tmgv0ijlyQsvvJB1KX+7tH3H/QAoLKEbgGZLXa61AEPbkqYFS6OU33HHHdk93Dt2Ma+trc22H3LIIdl+ABSegdQAANqxNA93mhZs2bJlceWVV8bvf//7eOONN7KfaT1tv+iii8zXDVAkWroBANq5NA/37Nmzs1HM06Bp9VILd9punm6A4hG6AQBKQArWJ554YjZKeRo0Ld3DnbqUa+EGKC6hGwCgRKSAbVowgLbFPd0AAACQE6EbAAAAciJ0AwAAQE6EbgAAAMiJ0A0AAAA5EboBAAAgJ0I3AAAA5EToBgAAgJwI3QAAAJAToRsAAAByInQDAABAToRuAAAAyEmnvJ64VNXV1WU/N23aVOyiAFCm6uug+jqpXKmTAWgP9bHQvYdee+217Ge/fv2KXRQAylyqkyorK6NcqZMBaA/1cYe6cr9Mvodqa2vjpZdeiv333z86dOhQ7OK0qas86Y+eNWvWxAEHHFDs4tCG+aywJ3xempaq7lTBH3rooVFRUb53iqmTm+b/Dc3ls0Jz+ay0rD7W0r2H0i/zsMMOK3Yx2qz0n9B/RJrDZ4U94fOys3Ju4a6nTt49/29oLp8VmstnZe/q4/K9PA4AAAA5E7oBAAAgJ0I3raJr164xa9as7Cfsjs8Ke8LnBfac/zc0l88KzeWz0jIGUgMAAICcaOkGAACAnAjdAAAAkBOhm1b1uc99Ls4444xiF4M26mc/+1kceeSR0bFjx/j7v//7XW57JwMHDox58+btdp80Z296bgrrQx/6ULPfxz3huwX2nP837Ir6uPSpj9sW83QDBfN3f/d3MWnSpPjiF78Y+++//y63Ub5efPHFOPzww+O3v/1tDB8+vGH79773vTAECUDrUB/zTtTHrUvoLkPbtm2LLl26FLsYlJnXX389XnnllRg3blwceuihu9wGTamsrCx2ESAX6mQKTX1MS6iP947u5WXSvWTq1KlZF5OePXtmX6hz586N97znPbHvvvtGv3794uKLL86+cOvdeuut0aNHj3jwwQfjmGOOif322y9OPfXUePnllxv22b59e0ybNi3b7+CDD46vfvWrO1352rp1a3bFtFevXtGtW7f4wAc+EL/+9a8bHl+6dGnW7Sid5/jjj4999tkn/uqv/ir74v/FL36RnfuAAw6Iz372s/HGG28U6DdW3u69997ss5Hei/S+jh07NjZv3pw9tmjRojj22GOz6SIOOeSQ7HNVb3efqfQ+1181T+9ves93tS35yU9+0nCe1HXtO9/5zm7L/Mc//jFOPvnk7DM2ZMiQePjhh3P7/fDO/vKXv2SfjVQxp++cq666quG7oaluhuk7JH3nJOmqepK+D9K+6furqe5saXv6bknfOwcddFD06dMnvva1rzV63nf6noNiUCfTXOpjWkp93HYI3WXihz/8YXYl/Ve/+lUsWLAgKioq4h//8R/j97//ffbYL3/5y+w/y45Shfrtb387br/99nj00Udj9erV8ZWvfKXh8fTFm/5jpi/+xx57LF599dW47777Gj1Hes70hZ3OsXLlyuxeofQHRtp3R+k/5w033BCPP/54rFmzJj7zmc9k9wjdeeed8fOf/zweeuihuP7663P+LZH+gDv77LPj85//fDz99NNZpfs3f/M32Rf0TTfdFJdccklMmTIlnnrqqbj//vuz97Pe7j5TJ5xwQqxatSr7d/o8pPPsatuKFSuy9/+ss87KzpM+G6mSqK8E3q62tjYrY/p8/+d//mf2+b788ssL8vuiaen979SpUyxfvjzrhpYq23/6p39q1rHpmOTf//3fs8/ET3/6092eJ1Xg6X2/7rrr4uqrr270B15zvuegGNTJvBP1Ma1BfdyGpHm6KW0f/OAH644//vjd7nPPPffUHXzwwQ3rP/jBD9JlsLpnn322Ydv8+fPrevfu3bB+yCGH1F133XUN62+99VbdYYcdVvfJT34yW3/99dfrOnfuXHfHHXc07LNt27a6Qw89tOG4//iP/8jO8+///u8N+8yZMyfb9txzzzVs+7u/+7u6cePGteC3QHOsWLEi+92/+OKLOz2W3rcZM2Y0+7ne/pn685//nD13es93t+2zn/1s3SmnnNLouS677LK6IUOGNKwPGDCg7rvf/W727wcffLCuU6dOdX/6058aHv/FL36RPe99993X7PLSet83xxxzTF1tbW3DtssvvzzbljT1vlRWVmbfOckLL7yQ7fPb3/620T4TJ05s+G6pP88HPvCBRvu8733vy87V3M8kFIM6meZQH9NS6uO2RUt3mRgxYkSj9XTV6iMf+Uj07ds361J03nnnxf/+7/826i7WvXv3OOKIIxrWU/el1MUsqampya56jR49uuHxdCVt5MiRDevPPfdcvPXWW3HiiSc2bOvcuXOMGjUqu2q7o+OOO67h3717987OPWjQoEbb6s9NfoYNG5Z9LlIXoDPPPDNuvvnm+POf/5z97l966aXssV1pzmeqOdJnY8fPTJLWU5e11H2yqf1TV6Ud70EbM2bMHp2T1vX+978/64q24/uxq/evJXb83nj7d1RrfiahtamTeSfqY1qD+rjtELrLROryseNohB//+Mez/yCpG1HqPjR//vyGAV12rIx3lP7T5jVa4Y7nSudp6typ2xL5StOEpO5A6d69dC9W6j541FFHxbp163Z7XHM/U9DU90gKAntjd98TPpO0Zepk3on6mLypjwtL6C5D6cOe/iOk+7/SFbDBgwdnV033RBqQIV3FSvdu7DhYQ3rueumKfP09azv+Z06DtqQKhLYpfVGmK9mzZ8/OpolI72Gq+NMAKkuWLMntM1UvDdSz42cmSevpOdMfIU3tn+453HFAoSeeeGKvzk3r2PF7of79ePe73529f+9617savVfpivuOV7rrR3Fu6VX41vxMQp7UyeyK+piWUh+3HaYMK0NpsI1U0aarpqeffnrDQC576tJLL41vfOMb2X/eo48+OhucYePGjY2u5F900UVx2WWXZaMZ9u/fPxtcIf2HPv/881v5VdFaX86pIv/oRz+ajW6b1tevX59VpGkAlQsvvDDb/td//dfx2muvZZ+dL3zhC632mUq+/OUvx/ve97645pprYvz48bFs2bJsQJ8bb7yxyf3TaK7pC3zixInxrW99KzZt2hQzZsxo4W+ClkgDPKVRlNOcr2mwpvS5qB/xNo2Mm97P1MUtVeRpkJ0dr5Cnz1caqXfx4sVx2GGHZSPg7s30JK35mYQ8qZNpivqY1qA+bju0dJfpfUKpMv7mN78ZQ4cOjTvuuCPmzJmzV1/G6Z6M9OWa/sOmezQ+9alPNdon/QHw6U9/Otvvve99bzz77LPZVCQHHnhgK74iWkuaCiaNinvaaadlFeeVV16ZfTmnSj29z2n02lTZpulDUlehdFW0NT9TSfqc/PjHP4677rore66ZM2dmo2CmKSqakkbETCP0vvnmm9m9iRdccEFUVVW16PdAy0yYMKHh/Ugj7KYwkEbZTdLnKd3zd9JJJ2XTDqXRl9P9ojveh5pGOP3+97+f3Rf4yU9+cq/K0JqfSciTOpmmqI9pDerjtqNDGk2t2IUAAACAUqSlGwAAAHIidAMAAEBOhG4AAADIidANAAAAORG6AQAAICdCNwAAAORE6AYAAICcCN0AAACQE6EbaOTFF1+MDh06xJNPPpmtL126NFvfuHFjtEVf+9rXYvjw4a3+vG39dQNQ2tTH0S5eNzSH0A0lZP369XHRRRdF//79o2vXrtGnT58YN25c/OpXv9rr5zzhhBPi5ZdfjsrKymz91ltvjR49erzjcWm/VEkec8wxOz12zz33ZI8NHDhwj8qSjvnZz362R8cAQKGpj4EddWq0BrRrn/70p2Pbtm3xwx/+MAYNGhTr1q2LJUuWxP/+7//u9XN26dIl+2Nhb+y7777xyiuvxLJly2LMmDEN22+55ZbsDxEAKEXqY2BHWrqhRKRuV//n//yf+OY3vxkf/vCHY8CAATFq1KiYPn16fOITn2h0dfqmm26Kv/7rv4599tkn+2Pg3nvvbVa3rvTvSZMmRU1NTbYtLak72a506tQpPvvZz8aiRYsatv3f//t/s+dJ29/uX/7lX+K9731vdOvWLSvX7Nmz4y9/+Uv2WP1V+E996lNNXpW//fbbs22pBeCss86K1157reGxrVu3xhe/+MXo1atX9twf+MAH4te//nWj4x944IEYPHhw9jtJv7/UrQ8A9pT6WH0Mbyd0Q4nYb7/9siV190qV2u5cddVV2VX4//qv/4pzzjknqxSffvrpZnVtmzdvXhxwwAFZF7e0fOUrX9ntMZ///Ofjxz/+cbzxxhsN3dxOPfXU6N27d6P90h8oEyZMiEsvvTT+8Ic/xPe///1s36qqquzx+kr5Bz/4QXbeHSvp5557Lnvd//Zv/5YtjzzySHzjG99oePyrX/1q/OQnP8laHFauXBlHHnlk1s3v1VdfzR5fs2ZN/M3f/E2cfvrp2b1zF1xwQVxxxRXv+PsAgLdTH6uPYSd1QMm499576w488MC6bt261Z1wwgl106dPr/uv//qvRvuk//YXXnhho22jR4+uu+iii7J/v/DCC9k+v/3tb7P1//iP/8jW//znP2frP/jBD+oqKyvfsSw77jd8+PC6H/7wh3W1tbV1RxxxRN2//Mu/1H33u9+tGzBgQMP+H/nIR+quvfbaRs9x++231x1yyCGNyn7fffc12mfWrFl13bt3r9u0aVPDtssuuyx7Tcnrr79e17lz57o77rij4fFt27bVHXrooXXXXXddtp5+T0OGDGn0vJdffnmj1w0AzaU+/v/Ux/D/aemGEpKulr/00ktx//33Z1evU7ex1D0sXaHe0Y73c9WvN+fK+t5KV9fTFfF0xXvz5s1x2mmn7bRPusp/9dVXN7QQpGXy5MnZVfT6q/K7krqx7b///g3rhxxySHbvWv1V97feeitOPPHEhsc7d+6cdfWrf83p5+jRo3f7OwKA5lIf/3/qY/j/hG4oMekeqVNOOSXrsvb444/H5z73uZg1a1ZRy5S6zD3xxBPZ/WbnnXdedm/Z273++uvZPWOpO1n98tRTT8Uf//jH7DXtTqq0d5TuMautrW311wEAzaU+Vh9DPaEbStyQIUOyq9k7ShXu29ebmkpkV6Onbt++fY/KcNBBB2WDx6Qr6+kqe1NSC8CqVauy+7vevlRUVDRU5nt67iOOOCIr847TtKQr7eketPS7SdJrX758+W5/RwDQEupj9THly5RhUCLSNCRnnnlmVoked9xxWfeu3/zmN3HdddfFJz/5yZ3m5Rw5cmQ2augdd9yRVXBp2pDmSF3H0lXwNPXJsGHDonv37tnyTlKXuhtvvDEOPvjgJh+fOXNmfPzjH8+mLvnbv/3brGJPXdx+97vfxde//vWGc6fzpq5pad7TAw88sFnTpKS5Ui+77LLsj430/Ol3krrInX/++dk+F154YXznO9/J9kmDtqxYsWKnLoAA0Bzq46apjylnWrqhRKR7rtJ9UN/97nfj5JNPjqFDh2Zd2tJ9WDfccEOjfVO3sbvuuiv7Y+C2226Lf/7nf264ytycEVNTpTh+/Ph417velVWYzZGm/thVBZ+k0UvTSKcPPfRQvO9974v3v//92WtJU63USxXxww8/HP369Yvjjz8+miuNnJrur0td6dIV/GeffTYefPDBhj8SUsWfRlNNI66mP1wWLFgQ1157bbOfHwDqqY93TX1MueqQRlMrdiGAwkn3V913331xxhlnFLsoAFC21MdQPrR0AwAAQE6EbgAAAMiJ7uUAAACQEy3dAAAAkBOhGwAAAHIidAMAAEBOhG4AAADIidANAAAAORG6AQAAICdCNwAAAORE6AYAAICcCN0AAAAQ+fh/o8qQaW47VgMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure, axes = plt.subplots(1,2,figsize=(10,5))\n",
    "ax_0 = sns.boxplot(x=\"method\",y=\"R2\",data=split_combo_df,ax=axes[0])\n",
    "ax_0.set_xlabel(\"Split Method\")\n",
    "ax_0.set_ylabel(\"$R^2$\")\n",
    "ax_1 = sns.boxplot(x=\"method\",y=\"MAE\",data=split_combo_df,ax=axes[1])\n",
    "ax_1.set_xlabel(\"Split Method\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59d445d7",
   "metadata": {},
   "source": [
    "[Tukey's Honest Significant Difference (HSD) test](https://en.wikipedia.org/wiki/Tukey%27s_range_test) allows us to evaluate whether there are statistically significant differences in model performance between the different data splitting methods. The table below summarizes the mean difference between each pair of methods (meandiff), the p-value adjusted for multiple comparisons (p-adj), the confidence interval boundaries, and whether the null hypothesis of equal means can be rejected (reject). The results indicate a signficant difference between the three splitting strategies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "8be966ac-9f1a-4e08-91ff-6dabbca683ff",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Multiple Comparison of Means - Tukey HSD, FWER=0.05  \n",
      "======================================================\n",
      "group1  group2  meandiff p-adj   lower   upper  reject\n",
      "------------------------------------------------------\n",
      "butina   random   0.4533    0.0  0.3457   0.561   True\n",
      "butina scaffold   0.3016    0.0  0.1939  0.4093   True\n",
      "random scaffold  -0.1517 0.0034 -0.2594 -0.0441   True\n",
      "------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "tukey_r2 = pairwise_tukeyhsd(endog=split_combo_df[\"R2\"], groups=split_combo_df[\"method\"], alpha=0.05)\n",
    "print(tukey_r2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a60efc29",
   "metadata": {},
   "source": [
    "We can use the same method to compare the MAE results.  In this case, all the differences are significant. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "30291155-6a0b-43b5-8419-afd47ec5caf0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Multiple Comparison of Means - Tukey HSD, FWER=0.05  \n",
      "======================================================\n",
      "group1  group2  meandiff p-adj   lower   upper  reject\n",
      "------------------------------------------------------\n",
      "butina   random  -0.1328    0.0 -0.1781 -0.0875   True\n",
      "butina scaffold  -0.0773 0.0003 -0.1226  -0.032   True\n",
      "random scaffold   0.0556 0.0123  0.0103  0.1009   True\n",
      "------------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "tukey_mae = pairwise_tukeyhsd(endog=split_combo_df[\"MAE\"], groups=split_combo_df[\"method\"], alpha=0.05)\n",
    "print(tukey_mae)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2b5f1bd",
   "metadata": {},
   "source": [
    "Define a small function to calculate an Analysis of Variance (ANOVA) on a dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "7ca22f17-f7bd-4a86-8880-c7dd0ef598c4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_anova(df_in, col, group_var=\"method\"):\n",
    "    res_list = []\n",
    "    for k,v in df_in.groupby(group_var):\n",
    "        res_list.append(v[col].values)\n",
    "    return f_oneway(*res_list)[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cfb8a6a",
   "metadata": {},
   "source": [
    "Find the splitting method that yielded the highest $R^2$ and lowest MAE."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f030fa04-57ed-4782-9d1c-63f2957b38ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('random', 'random')"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "best_r2 = split_combo_df.groupby(\"method\").mean().reset_index().sort_values(\"R2\",ascending=False).method.values[0]\n",
    "best_mae = split_combo_df.groupby(\"method\").mean().reset_index().sort_values(\"MAE\",ascending=True).method.values[0]\n",
    "best_r2,best_mae"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b77785de",
   "metadata": {},
   "source": [
    "The plot below is one of my favorites. It shows the mean value for each splitting method along with the associated confidence interval.  The colors in the plot designate priority and statistical significance.  \n",
    "- The best method is shown in blue. \n",
    "- When a value is shown in grey, we **cannot** reject the null hypothesis that the mean is equivalent to the value in blue.\n",
    "- When a value is shown in red, we **can** reject the null hypothesis that the mean is equivalent to the value in blue. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5e9c53cc-a456-4649-92ec-30f6582e30e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'ANOVA p=6.27e-09')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1wAAAIoCAYAAAB9IeqnAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAR0tJREFUeJzt3QmcHFWdOPDXmRwEwgwQroSEcIeE+5A73JeAghcefzmU9UAQWGVRl0VQwXALcoi6csgiLLd+1A2CHBJuuURB5BaQS46ZRCAhk/p/fuVOb08ySXeSrsxMzff7+XRSXV3dXe/N6/fqV+/Vq0qWZVkCAACg6QY1/yMBAAAIAi4AAICCCLgAAAAKIuACAAAoiIALAACgIAIuAACAggi4AAAACiLgAgAAKIiACwAAoCACLgAAgIIIuOh3zj///FSpVNKWW245z23i9XicccYZc7128cUX56/9/ve/n+u1O+64I33oQx9KK620Uho2bFhabbXV0he+8IX017/+tbrNq6++mgYPHpw+/elPz/P7p02bloYPH54+/OEPL/C+95YnnngifeITn0hjxoxJSy65ZFp33XXTt7/97fT222/Xfe+1116bPv7xj6c11lgjf+/48ePTV7/61fTWW2/Nte2//uu/pk033TQtt9xy+bYTJkxIJ5xwQpo+fXrT03TSSSelD37wg/nfM/I9vqcnsb6rzNQ+llhiiabvE1BO2qbiPPDAA3ld3tVurL/++un73/9+09qmW2+9tcc2oOsRbUkzzZgxI33ta19Lo0ePzv8eke833njjXNu999576Vvf+la+//F3j/9PPPHENGvWrKbuD8UbvBi+A5rqsssuyxube++9Nz355JNprbXWmue2p512Wjr00EPziraec845Jx155JF5hfblL385jRo1Kj322GPpP//zP9N///d/p1//+tdpm222SSuuuGLabbfd0s9//vM8GOnps6OSf/fdd+dq+BZk3xen559/Pm2xxRapra0tHX744Xmjdtddd6Xjjz8+3X///Xla5+fzn/983nBEelddddX0yCOPpHPPPTfPs2goo0Hpct9996VJkyalz3zmM3lA8+CDD6aTTz453XTTTel3v/tdGjSoeeeB/uM//iOtvPLKaZNNNkk33HBD3e1/8IMfpBEjRlSft7S0NG1fgHLTNhXjN7/5TfrABz6Q1+PHHXdcXkc/9dRT6YUXXqj73kbbpjjxd+mll871/lgX37/77rs3NU0HH3xwuvrqq9NRRx2V1l577TzY3muvvdItt9yStttuu+p2sd9XXXVV+uxnP5s233zzdPfdd+d5EIH2j370o6buEwXLoB95+umnsyi21157bbbCCitkJ5xwQo/bxTYbb7xx/v8ZZ5zR7bWLLrooX3/fffdV102dOjUbNGhQNmnSpOwf//hHt+2ffPLJbKWVVspGjRqVvfHGG/m6Sy+9NP+Myy+/vMfv33333bO2trbs3XffXeB97w0nnXRSvm9//OMfu60/8MAD8/Vd6Z6XW265Za51l1xySf7eH//4x3W///TTT8+3veuuu7JmeuaZZ/L/X3vttfzzjz/++B63i/XxemwHsKC0TcVob2/P0/ihD30o6+zsXOD3L2rbtNZaa2Vrr7121kz33HNP/v2nnXZadd0777yTrbnmmtnWW29dXXfvvffm2x133HHd3v/Vr341q1Qq2cMPP9zU/aJYhhTSr8RZuGWXXTbtvffe6aMf/Wj+fF623XbbtPPOO6dTTz01vfPOO/P93O985zv5sIFLLrlkrrOCa665Zv4ZL730UvrhD3+Yr4uhHUsttVT62c9+NtdnxbCO3/72t/n+xRCAhdn3OcWZx3322Sc/07bxxhvnPUMTJ07Mz1Y2Q0dHR/5/DFepFWdSo8dp6NCh833/jjvuONe6yKMQZ2IbSV+Yc5jHiy++mJ/Z6xpGs95666ULL7ywgRR1/9xGxfFQ5MU/j4sAGqNtKqZtinS88sor+ZC+aIv+8Y9/pNmzZzf8/kVpm7p6+/7f//t/c7325z//Oc+rGA0SaY7ep1/84hcN7VP0bMXoieh96xKfccghh+QjS2LESbj99tvz/2Oof614Hm1U9G7Sfwi46FeiIYix5xEAfPKTn8yvO4ohavMS1+ZEZR1DxeYlhl5EIxTD3FZfffUet4kx4NFA/fKXv8yfR4O277775sPU3njjjW7bRiXY2dk5VyW9oPs+p9g+9uP9739/mjx5cj5W/2Mf+1i3cd/REP39739v6BFjw+dslKLCf+ihh/IKP9IR+XbEEUfk6V1QL7/8cv7/8ssvP9drMf489uFvf/tb3lDH0L+ll146H9bYJf5uW221VT7UMIY5nn322fkwl9jHs846KxUhhuzEsMrYlxjKEfsAUI+2qZi2Ker/1tbW/ORbXH8VwwnjeQzHjKGRC2N+bdOc+RLmzK8//elPedsUAdvXv/71/Hq8yPf99tsvXXfddXW/P4bRr7POOnk6anW1f9EGd13nFWqH5IeuwDuG+9OPFNyDBk3z+9//Pu9ev/HGG/Pns2fPzsaMGZMdeeSRc20b2x122GH58k477ZStvPLK2dtvv93jsI2HHnoof97T59TacMMNs+WWW676/Fe/+lX+vh/+8Ifdtttqq62yVVZZpdvwhwXZ956MGzcuf/8111zTbahFDCXZZJNNug2hi+0aecw51OI73/lONnz48G7bHHvssdnCOuSQQ7KWlpbsL3/5y1yvxdDB2u8ZP378XPsT74/0/f3vf++2/hOf+EQ+JKbr79mIekMKzzrrrOzwww/PLrvssuzqq6/O/y6DBw/Oh5JEPgPMi7apuLYp0rbkkkvmjy9/+cv598T/sV20Bc1um7rMmjUrH8q4xRZbzPXaLrvskm2wwQbdhmVGvm2zzTYNDT9cb731sp133nmu9X/605/ydF1wwQX580hrPI9horXi9Vi//vrr1/0u+g6TZtBvxNmmGFq200475c9jmEWcVfuv//qv/AzTvCY4iDOJO+ywQ7rgggvyGfJ6mrUpRK/G/MTrXUPvQlxEu8IKK+RDHrqGBjzzzDP5Ra1HH310t8kfFnbfa8WFv11DIUKcHTvwwAPTKaeckp+xi8kh4tHTTEc92WijjeYaGrL99tunj3zkI2nkyJHpV7/6Vfrud7+bf2b0MC2IyJOf/OQn6ZhjjskvCJ5TDDmJ/YzhIXfeeWd+FrN2lsI4LrnmmmvS/vvvny/HWc8ue+yxR7riiivyC55jaE4zxAXptSIP4mxjnNmM2bviLCZAT7RNxbVN0S5ET98Xv/jF6qyE0Rs3c+bMfBhlzKTbUxuzsG1Tl+hZjB7If//3f++2PnoNb7755vx74+/T9TfqaptioqnojVtllVXm+dkxjLR2SGeXrllxu4aZxiQa48aNy/9m0au12WabpXvuuScde+yxeS9iveGo9DG9HfFBI+JsU5wxizNaTzzxRPVx5ZVX5md6brjhhnmeRZzzTOKinEVcdtllu62LXpG4oPmFF17In5944on5Zz344IMLve/zOou4/fbbz7X+Jz/5SVMmm4gLrKN36/nnn++2/uCDD87PLM7ZyzQ/v/vd77Illlgi22OPPbL33nuvofdEz1LkY/wtwiuvvFL3LGhc4B1eeumlbo+eer7q9XDNS5SZOJsJ0BNtU7FtU/QGxefcdttt3dbH81gfE2AU0TbFhFHRC/byyy/3OOHF/B4PPPBAnrdztk0zZsxYoB6uEBNZTZw4sfrZw4YNy84+++xsxRVXzDbaaKOG007v08NFvxBnlOLC4OjZiMec4izd/KZtjbNOcZ1SnBFbZpllur0W1wXF2aI//OEP83x/jKV+/PHH8wtja8V1PjHF7OWXX56fhYr/o/cmLh5u1r4viBif/9prrzW0bVzs2zUZRvTixJS7cQ+uWnHfk5iuNsac77rrrnU/8+GHH87fE/dIiQuDI18bEWcsDzjggDx/4uxm10XRkb8HHXRQj+/ZcMMNqxN71LrooovyKXebYezYsXNdBwHQRdtUbNsUvWdxzdScEzrFFPjhzTffbOgzF6Rtip6juBYr2rw5v7erbYo8jR6tnsTfLa6DnvO6u5jyPf7W0WZFL9ic4m/RleYuMVHUH//4x/Too4/maY2/YVzTFT2i0TtK/yHgol+Iij8q2PPOO2+u12I2pKgcY1jGnBeXdomKKSq6GOLwzW9+s9trcbFrDKeIxue5557Lu/DndOWVV+YNW8zGVCtuVhgzRcUwhbj/STQMc94gcVH3vUvMlhQnSGPIR5e//OUv3Wbj66mSn5euyj/E0ImYpWpOXRcvN3KTxbgvyp577pmnNe5xUns/q3oib6Mha29vz5/HcJgYJhONdL1Ab85hKtFANUPk9bPPPpsHogA90TYV2zbFMLqo47smzegSEy51tRXNbptitsEYKtjT7IQxsVIYMmTIfNumeH3OtqlrqGQEvZHGGAZaO3FGDBfser1W5GttuxZpiPaykZOg9CG93cUG9cRQi6WXXjr77Gc/2+Prd9xxR97VfsUVV8xz2Ea49dZbu90DpfZeJzE8IYZf7LjjjnMNSYt7lMSQj9p7ndT65je/mX/evvvum98bo+veTwu77wt6YXKkp/ZeHnHxcyOP2rTss88+2dChQ7PHH3+82/fut99+eb68+OKL1XXPPfdc9thjj3XbLoZLrLHGGtno0aO7pX9Ob775ZjZz5sx53ocrhqHUDmeMfXrkkUfm2v7VV1/NFkS9IYU9fd55552Xv+fMM89coO8CBgZtU/FtUwzPi8//1Kc+1e17P/nJT+YTGzWrbar1wQ9+MB9KP23atB5fj79FTFLyt7/9baHaprvvvnuu+3DFBBxxz68tt9xyvu+Nv9umm26a529HR0dD6aFvEHDR50WFH5XT9ddf3+PrMeNS3KzxAx/4wHwbtbDDDjtUx0LXNmohDqxjfVR6Md49Dv7/7d/+LVtmmWWy1tbW/AaUPYmZjro+c9ttt13kfZ9Xo7bOOuvk+/L1r389+973vpfPkhQN8ZQpU7JFFY16jFePceHf/va382Dj/e9/f77v//Iv/9JjHtaKseSx7phjjslnVKp9/OY3v6lud91112Vjx47N/vVf/zU7//zz89kBP/KRj+QHA5tvvnl1jHuIsfOR7mj44hqGmHFr8uTJ2cc+9rG5rleYl5/+9Kf57Ivf+MY38v2L6yXieTyeffbZ6nZx/VoEeHEj0kh7NOaxT3HAMOfNRgGCtqn4tilEUBj7uv/+++f1c7QB8Tzq9Wa1TV1ef/31bMiQIfOdATGutYo2aOTIkXmaf/SjH+Vtyl577ZVfT9eISEMEjPF3jLYtZjiM53NeqxbbdbV/EaBNmDAhv47rpptuauh76DsEXPR5UeHHha7zO/CNg+WoJLsmd5hXoxbTzc6rUeu6qDbOBi6//PL556266qrZ5z73uW4H5z153/vel39mBBGLuu/zatT23nvv/CLmqNCjwl133XWzq666KmuWuBg4gqw4Yxr7E43oSSedNNfFxT01avO7gDi27/Lkk0/mFyPHGccIciJv4gLi6HmaPn36XPsUk2fE3zGCtNinrkksooFrRO1BzPymHo6gMi5MjjO+8T1xYPO1r33NGURgnrRNi6dtilERJ5xwQv5dXfVzBHZzWpS2ac4p13/xi1/Md5+eeuqpvC3rai9juv0YKRK3FWlE9PgdffTR+fsjz+Lv1FOAesopp+T5GX+rCPKi96124hP6j0r809vDGoH5i3HwcbFv180tAaC3aZugMf93MwYAAACaSsAFAABQEAEXAABAQVzDBQAAUBA9XAAAAAURcAEAABRkcFEf3J/NmDEjf3SZPXt2euONN9LIkSNTpVLp1X0DGGhi5Pu0adPS6NGj06BBA/c8obYJoH+2SwKuHkyePDl961vf6u3dAKDG888/n8aMGZMGKm0TQP9sl0ya0cBZxPb29rTqqqvmmdra2tqr+wb0TdHbEHVEGDt2bOl6YnozfR0dHfl3vvXWW6mtrS0NVL3dNpW9jDeTvIJy61jAdkkPVw+GDRuWP+YUDZqAC+jJrFmz0nPPPZcvT5gwIQ0eXK7qtS+kb6APm+vttqkvlIH+Ql7BwFBpsF1yygUAAKAgAi4AAICCCLgAAAAKIuACAAAoiIALAACgIAIuAACAgpinFKAJ4j4773vf+6rLZVP29FGfMtA4eQXUEnABNEEcVK200kqprMqePupTBhonr4BaTrsAAAAURA8XQBPMnj07vfjii/nyKqusUrphRGVPH/UpA42TV0AtARdAkw6wHn744Xx51KhRpTvAKnv6qE8ZaJy8AmqpAQAAAAoi4AIAACiIgAsAAKAgAi4AAICCCLgAAAAKIuACAAAoiGnhAZogpn3edNNNq8tlU/b0UZ8y0Dh5BdQScAE0QRxUjR49OpVV2dNHfcpA4+QVUMtpFwAAgILo4QJogtmzZ6eXX345X1555ZVLN4yo7OmjPmWgcfIKqKUGAGjSAdYDDzyQP2K5bMqePupTBhonr4BaAi4AAICCCLgAAAAKIuACAAAoiIALAACgIAIuAACAggi4AAAACuI+XABNEPfZ2WijjarLZVP29FGfMtA4eQXUEnABNEEcVI0dOzaVVdnTR33KQOPkFVDLaRcAAICC6OECaILZs2en1157LV9eYYUVSjeMqOzpoz5loHHyCqilBgBo0gHWfffdlz9iuWzKnj7qUwYaJ6+AWgIuAACAggi4AAAACiLgAgAAKIiACwAAoCACLgAAgIIIuAAAAAriPlwATRD32Vl//fWry2VT9vRRnzLQOHkF1BJwATRBHFStttpqqazKnj7qUwYaJ6+AWk67AAAAFEQPF0ATZFmWXn/99Xx55MiRqVKppDIpe/qoTxlonLwCaunhAmiCzs7OdPfdd+ePWC6bsqeP+pSBxskroJaACwAAoCACLgAAgIIIuAAAAAoi4AIAACiIgAsAAKAgAi4AAICCuA8XQBMMGjQoTZgwobpcNmVPH/UpA42TV0AtARdAE8RB1ZprrpnKquzpoz5loHHyCqjltAsAAEBB9HABNEGWZam9vT1fbmtrS5VKJZVJ2dNHfcpA4+QVUEsPF0ATdHZ2pqlTp+aPWC6bsqeP+pSBxskroJaACwAAoCACLgAAgIIIuAAAAAoi4AIAACiIgAsAAKAgAi4AAICCuA8XQBMMGjQorb322tXlsil7+qhPGWicvAJqCbgAmiAOqsaPH5/Kquzpoz5loHHyCqjltAsAAEBB9HABNEGWZWn69On58ogRI1KlUkllUvb0UZ8y0Dh5BdTSwwXQBJ2dnem2227LH7FcNmVPH/UpA42TV0AtARcAAEBBBFwAAAAFEXABAAAURMAFAABQEAEXAABAQQRcAAAABXEfLoAmGDRoUFpjjTWqy2VT9vRRnzLQOHkF1BJwATRBHFRNnDgxlVXZ00d9ykDj5BVQy2kXAACAgujhAmiCLMvSO++8ky8PHz48VSqVVCZlTx/1KQONk1dALT1cAE3Q2dmZbr755vwRy2VT9vRRnzLQOHkF1BJwAQAAFETABQAAUBABFwAAQEEEXAAAAAURcAEAABREwAUAAFAQ9+ECaIK4z864ceOqy2VT9vRRnzLQOHkF1BJwATRBS0tL2mCDDVJZlT191KcMNE5eAbUMKQQAACiIHi6AJsiyLM2cOTNfHjp0aOmGEZU9fdSnDDROXgG19HABNEFnZ2e68cYb80csl03Z00d9ykDj5BVQS8AFAABQEAEXAABAQQRcAAAABRFwAQAAFETABQAAUBABFwAAQEHchwugCeI+O2PGjKkul03Z00d9ykDj5BVQS8AF0AQtLS1p4403TmVV9vRRnzLQOHkF1DKkEAAAoCB6uACaIMuy1NnZWT27XbZhRGVPH/UpA42TV0AtPVwATRAHV1OmTMkfXQdaZVL29FGfMtA4eQXUEnABAAAURMAFAABQEAEXAABAQQRcAAAABRFwAQAAFETABQAAUBD34QJogrjPzqhRo6rLZVP29FGfMtA4eQXUEnABNEHc3HSzzTZLZVX29FGfMtA4eQX02SGF119/fVprrbXyiuqoo46a57p6VltttXTWWWfNd5s44xSfDQAAMCB6uL7whS+kz3zmM+mII45ISy+99DzXASXR2ZnS7ben9NJLKcXwm0mT4tRwb+8VAD1RZ0P/DrimT5+eXn311bTHHnuk0aNHz3MdUBLXXpvSkUem9MIL/7duzJiUzj47pQ9/OPU3s2bNSlOmTMmX99xzzzR4cJ+pXpui7OmjPmVggOdVyeps6NNDCq+++uq0wQYbpOHDh6eRI0emXXfdNf3jH//IX7vwwgvTeuutl4YNG5ZfLHr44YdX33fmmWfm71tqqaXS2LFj05e+9KU8oAq33nprtfdq5513zof7zWtduOaaa6rfE8MHzzjjjPnu8xNPPJG23377tMQSS6SJEyemG2+8cUGTDTS74f7oR7s33OHFF/+5Pl4HoG9QZ8MiWaBTLi+99FL65Cc/mU499dT0oQ99KE2bNi3dfvvtKcuy9IMf/CB95StfSSeffHJ6//vfn9rb29Mdd9xRfe+gQYPS97///bT66qunp59+Og+4jjnmmHT++eenbbbZJj3++ONp/PjxeTAVz5dbbrke191///1p//33TyeccEL6+Mc/nu688878syL4O/jgg+fa59mzZ6cPf/jDaaWVVkr33HNPvl+NXgsGA9b/nkQpbEjKEUeklGVzvxbrYkavOIu6667FDVVZaqliPhegP9XH/aXObpS6nbIEXNFNHgHMuHHj8nXRaxVOPPHE9NWvfjUdGT+6//W+972vulwb5ESvVGz/xS9+MQ+4hg4dmlZcccX8tQiqVl555Xy5p3XRU7bLLruk4447Ln++zjrrpEcffTSddtppPQZcN910U/rzn/+cbrjhhuqwxO9+97t5UDgvM2bMyB9dOjo6FiSboP8bMaL3vjsa8DiL2tZW7HdAP6NtGqB6sz7uK3X2guwL9PchhRtttFEe7ESQ9bGPfSz9+Mc/Tm+++WZ+ndXf/va3/LV5icAnXl9llVXyoYIHHHBAev3119Pbb7+9QDv82GOPpW233bbbungewwY74yxMD9vHEMbaa8C23nrr+X7H5MmTU1tbW/UR7weA3qRtAhgAPVwxNXtc/xTD+H7zm9+kc845Jx177LHpt7/97Xzf9+yzz6Z99tknHXrooemkk07Ke6ymTp2aDjnkkDRz5sy05JJLpr7kG9/4Rj48svYsooaNAeV/r68sxO9+l9Jee9Xf7te/Tmn77YvbD+hntE0DVJH1cSPU2bDIFnjanJi8InqU4vHNb34zH1oYQVgME4zAa6eddprrPXHdVVxLFZNbxLVc4corr1yoHZ4wYUK3a8NCPI+hhREQ9rT9888/nw+H7Lrr+9133z3f74jJOOIBA1aR4+B33/2fM1vFxdY9Df+I6wHi9diut68HgD5E2zRA9fZ1SepsWLwBV0w6EUHV7rvvnl9fFc9fe+21PKiJSSzimqxYH9dHxYQaEQh9+ctfzm9c/N577+U9Yh/4wAfy9RdccMFC7XBcJxbXhn3nO9/JJ82466670rnnnptfC9aTmEUxgrGDDjoov84rzghGrxzQS6JBjmmEY2araKhrG/B4HuLG5f2s4Y6TUV3XncZy2ZQ9fdSnDAzQvCppnQ2LVbYAHn300WyPPfbIVlhhhWzYsGHZOuusk51zzjnV1y+44IJs/Pjx2ZAhQ7JRo0ZlX/7yl6uvnXnmmfm64cOH55/x05/+NH6x2Ztvvpm/Hv/H81tuuaX6np7WhauvvjqbOHFi/j2rrrpqdtppp3V7fdy4cdn3vve96vPHH38822677bKhQ4fm+zxlypT8c6+77rqG0t3e3p5vH/8DTXLNNVk2Zkw03f/3GDv2n+uhhjq4Z/KFxUqdDQtd/1bin8Ub4vU/0SsWFyjHlPKtra29vTtQHjHRze23xxSoKcWQ30mTnCVlLurgnskXFjt1NixU/VuCW58D/VY01Dvu2Nt7AUAj1NmwUARcAE0Q9yiMCYTCbrvtlgYPLlf1Wvb0UZ8y0Dh5BdRSAwA0SU/3AiyTsqeP+pSBxskrYKFufAwAAEDjBFwAAAAFEXABAAAURMAFAABQEAEXAABAQcxSCNAElUolLbfcctXlsil7+qhPGWicvAJqVbIsy7qtYZHvJg1A86iDeyZfAPpH/WtIIQAAQEEEXAAAAAVxDRdAE8yaNSvdfPPN+fLOO++cBg8uV/Va9vRRnzLQOHkF1FIDADTJzJkzU5mVPX3Upww0Tl4BXQwpBAAAKIiACwAAoCACLgAAgIIIuAAAAAoi4AIAACiIWQoBmqBSqeR3ne9aLpuyp4/6lIHGySugViXLsqzbGubS0dGRV5zt7e2ptbW1t3cHYEBRB/dMvgD0j/rXkEIAAICCCLgAAAAK4hougCbo7OxMt956a7684447ppaWllQmZU8f9SkDjZNXQC0BF0ATxOWw77zzTnW5bMqePupTBhonr4BahhQCAAAURMAFAABQEAEXAABAQQRcAAAABRFwAQAAFMQshQBNUKlU0ogRI6rLZVP29FGfMtA4eQXUqmTmK62ro6MjtbW1pfb29tTa2trbuwMwoKiDeyZfAPpH/WtIIQAAQEEEXAAAAAVxDRdAE3R2dqbbb789X540aVJqaWlJZVL29FGfMtA4eQXUEnABNEFcDjt9+vTqctmUPX3Upww0Tl4BtQwpBAAAKIiACwAAoCACLgAAgIIIuAAAAAoi4AIAACiIWQoBmqBSqaThw4dXl8um7OmjPmWgcfIKqFXJzFdaV0dHR2pra0vt7e2ptbW1t3cHYEBRB/dMvgD0j/rXkEIAAICCCLgAAAAK4hougCbo7OxMd955Z768zTbbpJaWllQmZU8f9SkDjZNXQC0BF0ATxOWwMZa7a7lsyp4+6lMGGievgFqGFAIAABREwAUAAFAQARcAAEBBBFwAAAAFEXABAAAUxCyFAE0ydOjQVGZlTx/1KQONk1dAl0pmvtK6Ojo6UltbWz7Fa2tra2/vDsCAog7umXwB6B/1ryGFAAAABRFwAQAAFMQ1XABN0NnZme655558ecstt0wtLS2pTMqePupTBhonr4BaAi6AJojLYd94443qctmUPX3Upww0Tl4BtQwpBAAAKIiACwAAoCACLgAAgIIIuAAAAAoi4AIAACiIWQoBmqTsUz+XPX3Upww0Tl4BXSqZ+Urr6ujoSG1tbam9vT21trb29u4ADCjq4J7JF4D+Uf8aUggAAFAQARcAAEBBXMMF0ASdnZ3p/vvvz5c322yz0l2/Ufb0UZ8y0Dh5BdQScAE0QVwO++qrr1aXy6bs6aM+ZaBx8gqoZUghAABAQQRcAAAABRFwAQAAFETABQAAUBABFwAAQEEEXAAAAAWpZOYrraujoyO1tbWl9vb21Nra2tu7AzCgqIN7Jl8A+kf9q4cLAACgIAIuAACAggwu6oMBBpLOzs700EMP5csbb7xxamlpSWVS9vRRnzLQOHkF1NLDBdAEcTnsSy+9lD/KeGls2dNHfcpA4+QVUEvABQAAUBABFwAAQEEEXAAAAAURcAEAABREwAUAAFAQARcAAEBBKpn5Suvq6OhIbW1tqb29PbW2tvb27gB9UFSlce+dEPfcqVQqqUx6M33q4L6RL2Uv480kr6DcOhaw/nXjY4AmiAOqwYPLW6WWPX3Upww0Tl4BtQwpBAAAKIjTLwBNEMOHHnnkkXx5gw02yIcRlUnZ00d9ykDj5BVQSw8XQJOu2XjhhRfyRxkvjS17+qhPGWicvAJqCbgAAAAKIuACAAAoiIALAACgIAIuAACAggi4AAAACiLgAgAAKEglM19pXR0dHamtrS21t7en1tbW3t4doA+KqnTmzJn58tChQ1OlUkll0pvpUwf3jXwpexlvJnkF5daxgPWvGx8DNEEcUA0bNiyVVdnTR33KQOPkFVDLkEIAAICC6OECaILOzs706KOP5ssTJ05MLS0tqUzKnj7qUwYaJ6+AWnq4AJp0zcZzzz2XP8p4aWzZ00d9ykDj5BVQS8AFAABQEAEXAADAQLyG6+CDD05vvfVWuv7663t7VwD6lM7OlG6/PaWXXkpp1KiUJk1KyWUiDER+C0Bf16cDLgDmdu21KR15ZEovvPB/68aMSenss1P68Id7c89g8fJbAAbEkMKuG/sBsHgOMD/60e4HmOHFF/+5Pl6HgcBvAShtwLXjjjumww8/PB111FFp+eWXT3vssUc688wz0wYbbJCWWmqpNHbs2PSlL30pTZ8+vfqeiy++OC2zzDLphhtuSBMmTEgjRoxIe+65Z3op+v9rplD9yle+km83cuTIdMwxx8w1s8+MGTPSEUcckVZcccW0xBJLpO222y7dd9991ddvvfXW/GaD8T2bbLJJGj58eNp5553Tq6++mv7nf/4n/+64G/SnPvWp9Pbbby98rgHMR1Qv//hH8x8dHSkdcUTMgDb3d3ati7P9sV2zv1uVyeIo4/3ht+D3AiyWIYWXXHJJOvTQQ9Mdd9yRP49g5vvf/35affXV09NPP50HXBEwnX/++dX3RIBz+umnp0svvTQNGjQoffrTn05HH310uuyyy/LXzzjjjDwwu/DCC/PAKJ5fd911ecDUJT7zmmuuyb9/3Lhx6dRTT80DvieffDItt9xy1e1OOOGEdO6556Yll1wy7b///vkj7vj+s5/9LA8EP/ShD6Vzzjknfe1rX+sxfRHYxaNLR9TYAPMR99mJ+mr11VN67bWWHg8EixbfGWf729qa/9mVSktaYYWd0zPP/DOtLH693Tb1hTLeF34LjfB7AbrJFtAOO+yQbbLJJvPd5qqrrspGjhxZfX7RRRdFtZw9+eST1XXnnXdettJKK1Wfjxo1Kjv11FOrz997771szJgx2b777ps/nz59ejZkyJDssssuq24zc+bMbPTo0dX33XLLLfn33HTTTdVtJk+enK976qmnquu+8IUvZHvsscc89//444/P3zPno729vaE8Agaufx7qlffRG6LuVQf3nbapt8tgf3oA5bSg7dJC9XBtttlm3Z7fdNNNafLkyenPf/5zfsZt1qxZ6d133817taKXKcT/a665ZvU9o0aNyof6hfb29nx44ZZbbll9ffDgwWnzzTevDit86qmn0nvvvZe23Xbb6jZDhgxJW2yxRXrssce67c+GG25YXV5ppZXy715jjTW6rbv33nvnmb5vfOMb+fDGLpGmGCoJUE/NaOqm+93vUtprr/rb/frXKW2/fXH7Qe/oK21TkWW8UX4LQH+yUAFXXKvV5dlnn0377LNPPsTwpJNOyof2TZ06NR1yyCH5hBpdAVcER7XiWqui7r5e+13xPT199+zZs+f5/hh+GA+ARkWdEiedwrrrrpsPnW623Xf/5wxsMSlAT9VnpfLP12O7Zo9iWhzpY/56u23qS2WgN38L/S2vgN63yDXA/fffn1cscc3VVlttldZZZ530t7/9bYE+o62tLe/xuueee6rropcsPrtL9I4NHTq0et1YiB6vmDRj4sSJi5oMgEUS9WBcwxqP+Z3QWRRx4BjTXXcdUNbqen7WWcUcYC6O9NG39aUy0Ju/hf6WV0AJAq611lorD3xiEoqoWGJSjAsuuGCBP+fII49MJ598cn6T4zgrFBNvxE2Pa3vVohft3/7t39KUKVPSo48+mj73uc/lwxajNw1gIIh7C119dUqrrNJ9fZzNj/XuPcRA4bcADJgbH2+00Ub5tPCnnHJKPr58++23z6/nOvDAAxfoc7761a/m13EddNBBedf7Zz/72Xw2wbi+q0sEZHGm6IADDkjTpk3Lr/GKKeCXXXbZRU0GQL8RB5L77pvS7benFHfXGDUqpUmTeu9sPvQWvwWgP6jEzBm9vRN9XVyYHMMeI/iL+3gBzCmGQUfve4j7DMbEP2XSm+lTB/eNfCl7GW8meQXl1rGA9a+rOAEAAAoi4AIAACiIgAsAAKAgBhUDNEFLS0vaYYcdqstlU/b0UZ8y0Dh5BdQScAE0QdxQfemll05lVfb0UZ8y0Dh5BdQypBAAAKAgergAmiDuEfjEE0/ky2uvvXZ+P8EyKXv6qE8ZaJy8AmoJuACafIC15pprlu4Aq+zpoz5loHHyCqilBgAAACiIgAsAAKAgAi4AAICCCLgAAAAKIuACAAAoiIALAACgIKaFB2iClpaWtN1221WXy6bs6aM+ZaBx8gqoJeACaIJKpZKWWWaZVFZlTx/1KQONk1dALUMKAQAACqKHC6AJZs+enZ555pl8efXVV0+DBpXrfFbZ00d9ykDj5BVQS8AF0KQDrMceeyxfHjduXOkOsMqePupTBhonr4BaagAAAICCCLgAAAAKIuACAAAoiIALAACgIAIuAACAggi4AAAACmJaeIAmaGlpSVtttVV1uWzKnj7qUwYaJ6+AWgIugCaoVCpp+eWXT2VV9vRRnzLQOHkF1DKkEAAAoCB6uACaYPbs2emvf/1rvrzqqqumQYPKdT6r7OmjPmWgcfIKqCXgAmjSAdYf//jHfHnMmDGlO8Aqe/qoTxlonLwCaqkBAAAACiLgAgAAKIiACwAAoCACLgAAgIIIuAAAAAoi4AIAACiIaeEBmiCmfX7f+95XXS6bsqeP+pSBxskroJaAC6AJ4qBqpZVWSmVV9vRRnzLQOHkF1HLaBQAAoCB6uACaYPbs2enFF1/Ml1dZZZXSDSMqe/qoTxlonLwCagm4AJp0gPXwww/ny6NGjSrdAVbZ00d9ykDj5BVQSw0AAABQEAEXAABAQQRcAAAABRFwAQAAFETABQAAUBABFwAAQEFMCw/QBDHt86abblpdLpuyp4/6lIHGySugloALoAnioGr06NGprMqePupTBhonr4BaTrsAAAAURA8XQBPMnj07vfzyy/nyyiuvXLphRGVPH/UpA42TV0AtNQBAkw6wHnjggfwRy2VT9vRRnzLQOHkF1BJwAQAAFETABQAAUBABFwAAQEEEXAAAAAURcAEAABREwAUAAFAQ9+ECaIK4z85GG21UXS6bsqeP+pSBxskroJaAC6AJ4qBq7NixqazKnj7qUwYaJ6+AWk67AAAAFEQPF0ATzJ49O7322mv58gorrFC6YURlTx/1KQONk1dALTUAQJMOsO677778EctlU/b0UZ8y0Dh5BdQScAEAABREwAUAAFAQARcAAEBBBFwAAAAFEXABAAAURMAFAABQEPfhAmiCuM/O+uuvX10um7Knj/qUgcbJK6CWgAugCeKgarXVVktlVfb0UZ8y0Dh5BdRy2gUAAKAgergAmiDLsvT666/nyyNHjkyVSiWVSdnTR33KQOPkFVBLDxdAE3R2dqa77747f8Ry2ZQ9fdSnDDROXgG1BFwAAAAFEXABAAAURMAFAABQEAEXAABAQQRcAAAABRFwAQAAFMR9uACaYNCgQWnChAnV5bIpe/qoTxlonLwCagm4AJogDqrWXHPNVFZlTx/1KQONk1dALaddAAAACqKHC6AJsixL7e3t+XJbW1uqVCqpTMqePupTBhonr4BaergAmqCzszNNnTo1f8Ry2ZQ9fdSnDDROXgG1BFwAAAAFEXABAAAURMAFAABQEAEXAABAQQRcAAAABRFwAQAAFMR9uACaYNCgQWnttdeuLpdN2dNHfcpA4+QVUEvABdAEcVA1fvz4VFZlTx/1KQONk1dALaddAAAACqKHC6AJsixL06dPz5dHjBiRKpVKKpOyp4/6lIHGySuglh4ugCbo7OxMt912W/6I5bIpe/qoTxlonLwCagm4AAAACiLgAgAAKIiACwAAoCACLgAAgIIIuAAAAAoi4AIAACiI+3ABNMGgQYPSGmusUV0um7Knj/qUgcbJK2CxBFw77rhj2njjjdNZZ53V1M89+OCD01tvvZWuv/76pn4uwKKIg6qJEyemsip7+qhPGWicvAL6RQ/Xs88+m1ZfffX04IMP5oFbl7PPPju/gztNFjdmvP32lF56KaVRo1KaNCmllpbe3isAAPoCx4rlC7jmpa2trbd3oXyuvTalI49M6YUX/m/dmDER3ab04Q/35p5BvxEngt555518efjw4alSqaQyKXv6qE8ZaJy8onQcKy6SQgcWz5o1Kx1++OF5kLT88sun4447rto7FZXPnMMCl1lmmXTxxRfny9G7FTbZZJN82xii2DWkcL/99qu+J9YfccQR6ZhjjknLLbdcWnnlldMJJ5zQ7XPPPPPMtMEGG6SllloqjR07Nn3pS19K06dPLzLp/esH9NGPdv8BhRdf/Of6eB2oq7OzM9188835I5bLpuzpoz5loHHyilJxrNi3e7guueSSdMghh6R77703/f73v0+f//zn06qrrpo+97nP1X1vvGeLLbZIN910U1pvvfXS0KFD5/s9X/nKV9I999yT7rrrrjwo23bbbdNuu+1WHUv9/e9/Pw/inn766TzgigDt/PPPT33eP/5R3GdHI3DEEXEqbu7XYl2ckYuzGbvuWlyX8VJLFfO5AAD9WZHHgP3pWLEEx5OFBlzRm/S9730v76EaP358euSRR/LnjQRcK6ywQv7/yJEj816r+dlwww3T8ccfny+vvfba6dxzz02//e1vqwHXUUcdVd12tdVWSyeeeGL64he/OM+Aa8aMGfmjS0dHR+o1I0b03nfHDynOZhQ5jNP1eAAN6VNtE1DuY8C+dKxYguPJQocUbrXVVt3GLW+99dbpiSeeaHr3egRctUaNGpVeffXV6vPoJdtll13SKquskpZeeul0wAEHpNdffz29/fbbPX7e5MmT82GQXY8IHAGgN2mbAPqnXps0IwKxOWcbfO+99xbqs4YMGTLXZ8+ePbs62+E+++yTDj300HTSSSfl13lNnTo1H+o4c+bMtOSSS871ed/4xjfyIYq1ZxF7rWEr8lqz3/0upb32qr/dr3+d0vbbF7cfANTVp9omoHh9Yb4Bx4p9P+CKa6pq3X333fmQv5aWlnzI4EsxreT/ip6v2h6nrmu2FrU37P7778+DrzPOOKN688Err7xyvu8ZNmxY/ugTihyTuvvu/5xhJi567KkrNnon4/XYzrSfAL2qT7VNwMC4LsmxYt8fUvjXv/41Pxv3+OOPp8svvzydc8456ci4sC6ltPPOO+fXWsV9tmJCjbimqranasUVV8ynUp0yZUp65ZVXUnt7+0Ltw1prrZX3nMV3x4QZl156abrgggualsZ+LX4YMZ1nmHPK2q7nceNqPyAAgIHHsWLfD7gOPPDA/D4UMdvgYYcdlgdbMVNhiB6nGAoxadKk9KlPfSodffTR3Yb3DR48OJ9Z8Ic//GEaPXp02nfffRdqHzbaaKN8WvhTTjklrb/++umyyy7Lx8Hzv+LeCVdfndIqq3RfH2crYr17K0BDYijzuHHj8kcZ77lT9vRRnzLQOHlFqThWXGSVbM4LqZhLjJOPC5Sjl621tTWVkruHA33UgKiDF4J8ARYrx4oLXf/22qQZ9DHxg/nfm0sDAEA3jhUXmoALoAlisEDMfNo16U/ZhhGVPX3Upww0Tl4Bi+0aLoCBImZUvfHGG/NHs+812BeUPX3Upww0Tl4BtQRcAAAABRFwAQAAFETABQAAUBABFwAAQEEEXAAAAAURcAEAABTEfbgAmiDuszNmzJjqctmUPX3Upww0Tl4BtQRcAE3Q0tKSNt5441RWZU8f9SkDjZNXQC1DCgEAAAqihwugCbIsS52dndWz22UbRlT29FGfMtA4eQXU0sMF0ARxcDVlypT80XWgVSZlTx/1KQONk1dALQEXAABAQQRcAAAABRFwAQAAFETABQAAUBABFwAAQEEEXAAAAAVxHy6AJoj77IwaNaq6XDZlTx/1KQONk1dALQEXQBPEzU0322yzVFZlTx/1KQONk1dALUMKAQAACiLgAgAAKIghhQBNMGvWrDRlypR8ec8990yDB5erei17+qhPGWicvAJq6eECAAAoiIALAACgIAIuAACAggi4AAAACiLgAgAAKIiACwAAoCDmKQVogkqlklZcccXqctmUPX3Upww0Tl4BtQRcAE3Q0tKStthii1RWZU8f9SkDjZNXQC1DCgEAAAoi4AIAACiIIYUATTBr1qx044035su77bZbGjy4XNVr2dNHfcpA4+QVUEsNANAknZ2dqczKnj7qUwYaJ6+ALoYUAgAAFETABQAAUBABFwAAQEEEXAAAAAURcAEAABTELIUATVCpVNJyyy1XXS6bsqeP+pSBxskroFYly7Ks2xrm0tHRkdra2lJ7e3tqbW3t7d0BGFDUwT2TLwD9o/41pBAAAKAgAi4AAICCuIYLoAlmzZqVbr755nx55513ToMHl6t6LXv6qE8ZaJy8AmqpAQCaZObMmanMyp4+6lMGGievgC6GFAIAABREwAUAAFAQARcAAEBBBFwAAAAFEXABAAAUxCyFAE1QqVTyu853LZdN2dNHfcpA4+QVUKuSZVnWbQ1z6ejoyCvO9vb21Nra2tu7AzCgqIN7Jl8A+kf9a0ghAABAQQRcAAAABXENF0ATdHZ2pltvvTVf3nHHHVNLS0sqk7Knj/qUgcbJK6CWgAugCeJy2Hfeeae6XDZlTx/1KQONk1dALUMKAQAACiLgAgAAKIiACwAAoCACLgAAgIIIuAAAAApilkKAJqhUKmnEiBHV5bIpe/qoTxlonLwCalUy85XW1dHRkdra2lJ7e3tqbW3t7d0BGFDUwT2TLwD9o/41pBAAAKAgAi4AAICCuIYLoAk6OzvT7bffni9PmjQptbS0pDIpe/qoTxlonLwCagm4AJogLoedPn16dblsyp4+6lMGGievgFqGFAIAABREwAUAAFAQARcAAEBBBFwAAAAFEXABAAAUxCyFAE1QqVTS8OHDq8tlU/b0UZ8y0Dh5BdSqZOYrraujoyO1tbWl9vb21Nra2tu7AzCgqIN7Jl8A+kf9a0ghAABAQQRcAAAABXENF0ATdHZ2pjvvvDNf3mabbVJLS0sqk7Knj/qUgcbJK6CWgAugCeJy2BjL3bVcNmVPH/UpA42TV0AtQwoBAAAKIuACAAAoiIALAACgIAIuAACAggi4AAAACmKWQoAmGTp0aCqzsqeP+pSBxskroEslM19pXR0dHamtrS2f4rW1tbW3dwdgQFEH90y+APSP+teQQgAAgIIIuAAAAAriGi6AJujs7Ez33HNPvrzlllumlpaWVCZlTx/1KQONk1dALQEXQBPE5bBvvPFGdblsyp4+6lMGGievgFqGFAIAABREwAUAAFAQARcAAEBBBFwAAAAFEXABAAAUxCyFAE1S9qmfy54+6lMGGievgC6VzHyldXV0dKS2trbU3t6eWltbe3t3AAYUdXDP5AtA/6h/DSkEAAAoiIALAACgIK7hAmiCzs7OdP/99+fLm222Wemu3yh7+qhPGWicvAJqCbgAmiAuh3311Very2VT9vRRnzLQOHkF1DKkEAAAoCACLgAAgIIIuAAAAAoi4AIAACiIgAsAAKAgZinswYwZM/JHl7iLdNddpQF6MmvWrPT2229X64rBg8tVvfZm+rrq3oE+21tvt01lL+PNJK+g3DoWsF1SA/Rg8uTJ6Vvf+tZc68eOHdsr+wNAStOmTUttbW1poNI2AfTPdqmSDfRThg2cRZw9e3Z644030siRI1OlUlmg6Dcawueffz61trYWtLeLV9nSJD19X9nSJD0LLpqpaNRGjx6dBg0auCPhF7RtKltZ603ysjnkY/PIy97NxwVtl/Rw9WDYsGH5o9Yyyyyz0J8Xf8Cy/RjKlibp6fvKlibpWTADuWdrUdumspW13iQvm0M+No+87L18XJB2aeCeKgQAACiYgAsAAKAgAq4CxdCP448/fq4hIP1Z2dIkPX1f2dIkPSwu/jbNIy+bQz42j7zsX/lo0gwAAICC6OECAAAoiIALAACgIAIuAACAggi4AAAACiLgWkTnnXdeWm211dISSyyRttxyy3TvvffOd/urrroqrbvuuvn2G2ywQfr1r3+d+nOa/vSnP6WPfOQj+faVSiWdddZZqT+n58c//nGaNGlSWnbZZfPHrrvuWvdv2pfTc+2116bNN988vznqUkstlTbeeON06aWXpv7+O+pyxRVX5OVuv/32S/01PRdffHGehtpHvK8//33eeuutdNhhh6VRo0blMz+ts846fbKu64+a/VuJebO++c1v5n+r4cOH53XeE088kcqu2fl48MEHz/U73nPPPdNA0Oz6TplsTj4O1DJ5XgHt1cLWF93ELIUsnCuuuCIbOnRoduGFF2Z/+tOfss997nPZMsssk73yyis9bn/HHXdkLS0t2amnnpo9+uij2X/8x39kQ4YMyR555JGsv6bp3nvvzY4++ujs8ssvz1ZeeeXse9/7XtaXLGh6PvWpT2XnnXde9uCDD2aPPfZYdvDBB2dtbW3ZCy+8kPXH9Nxyyy3Ztddem5e3J598MjvrrLPyMjhlypSsr1jQNHV55plnslVWWSWbNGlStu+++2b9NT0XXXRR1tramr300kvVx8svv5z11/TMmDEj23zzzbO99tormzp1av53uvXWW7OHHnpose972RTxWzn55JPzOu7666/PHn744eyDH/xgtvrqq2fvvPNOVlZF5ONBBx2U7bnnnt1+x2+88UZWdkXUd8pkc/JxIJbJKwporxa2vpiTgGsRbLHFFtlhhx1Wfd7Z2ZmNHj06mzx5co/b77///tnee+/dbd2WW26ZfeELX8j6a5pqjRs3rs8FXIuSnjBr1qxs6aWXzi655JKsDOkJm2yySR7s9xULk6b4u2yzzTbZf/7nf+aNSl8KuBY0PdFwxsFFX7Wg6fnBD36QrbHGGtnMmTMX414ODM3+rcyePTs/UXbaaadV17311lvZsGHD8pNoZVVEndPX6qH+Wt8pk81rNwZimdyigPaqGcddwZDChTRz5sx0//33513dXQYNGpQ/v+uuu3p8T6yv3T7sscce89y+P6SpL2tGet5+++303nvvpeWWWy719/TECZbf/va36fHHH0/bb7996gsWNk3f/va304orrpgOOeSQ1JcsbHqmT5+exo0bl8aOHZv23XfffKhuf03PL37xi7T11lvnQzRWWmmltP7666fvfve7qbOzczHuefkU8Vt55pln0ssvv9ztM9va2vIhM/2xzu/tOufWW2/Ntxk/fnw69NBD0+uvv57KrIj6TplsbrsxkMrkzALaq2YeFwu4FtLf//73/A8Sf6Ba8Twqi57E+gXZvj+kqS9rRnq+9rWvpdGjR88VKPen9LS3t6cRI0akoUOHpr333judc845abfddkt9wcKkaerUqeknP/lJfr1dX7Mw6YmG8MILL0w///nP03/913+l2bNnp2222Sa98MILqT+m5+mnn05XX311/r4YB3/cccelM844I5144omLaa/LqYjfStf7ylLn92adE9fG/PSnP81Pap1yyinptttuS+9///tLfaKhiPpOmWxeuzHQyuTfC2ivmnlcPHiBUwQDxMknn5xfIB1niPraJAYLYumll04PPfRQfjYsKt6vfOUraY011kg77rhj6m+mTZuWDjjggPzAZ/nll09lEGfX4tElGs0JEyakH/7wh+k73/lO6m+i4Y8zqj/60Y9SS0tL2myzzdKLL76YTjvttHT88cf39u4NGGX8rfTlfPzEJz5RXY4JsTbccMO05ppr5u3HLrvsspj2tu8rW33Xl/NRmexb7ZWAayFFxRt/nFdeeaXb+ni+8sor9/ieWL8g2/eHNPVli5Ke008/PQ+4brrppryS6s/pie7vtdZaK1+OWQofe+yxNHny5D4RcC1omp566qn07LPPpg984APdKswwePDgfLhkNCj9+Tc0ZMiQtMkmm6Qnn3wy9baFSU/M9BRpiPd1iQOBOBsYwzOip5W+8Vvpel98Rvzdaj8z6ooyWlx1TpzUiu+K33FZD26LqO+UyeLajbKXyeULaK+aeVxsSOFCioOGiISjx6C2Eo7ntWcdasX62u3DjTfeOM/t+0Oa+rKFTc+pp56anyGaMmVKPqV6X9Gsv0+8Z8aMGak/piluqfDII4/kPXZdjw9+8INpp512ypdjLHt//xvF8IVIY+3BRn9Kz7bbbps36F0HpeEvf/lLnh7BVt/6ray++ur5QUPtZ3Z0dKR77rmnX9b5fanOiaFdcb1MX/gd96f6Tpksrt0oe5kcWkB71dTj4gWaYoNuYqrImDnn4osvzqfd/vznP59PFdk1NecBBxyQff3rX+82LfzgwYOz008/PZ9y/Pjjj++T08IvSJpiSs2YQj0eo0aNyqeIj+Unnngi64/pieloY/rPq6++uttUqtOmTcv6Y3q++93vZr/5zW+yp556Kt8+yl6UwR//+MdZX7GgaerrMzEtaHq+9a1vZTfccEP+N7r//vuzT3ziE9kSSyyRTz/bH9Pz17/+NZ/Z8/DDD88ef/zx7Je//GW24oorZieeeGIvpqIcivitRJ0Xn/Hzn/88+8Mf/pC/PhCm4G5mPkb7EG3fXXfdlU8rfdNNN2Wbbrpptvbaa2fvvvtuVmZF1HfK5KLn40Atk1cU0F7V+8xGCbgW0TnnnJOtuuqq+UF6TB159913V1/bYYcd8oq51pVXXpmts846+fbrrbde9qtf/Srrz2mKH3LE7XM+Yrv+mJ6Y2r6n9ERw3B/Tc+yxx2ZrrbVWXhEvu+yy2dZbb51XHv39d9SXA64FTc9RRx1V3XallVbK7wfywAMPZP3573PnnXfmt7yIRiqm3D3ppJPyabXpe7+VmIb7uOOOy8te/L122WWX/MCj7JqZj2+//Xa2++67ZyussEJ+EjXakbhXT1+6n15/qu+UyUXPx4FcJs8poL2a32c2qhL/LEzXHQAAAPPnGi4AAICCCLgAAAAKIuACAAAoiIALAACgIAIuAACAggi4AAAACiLgAgAAKIiACwAAoCACLgAAgIIIuGCAef7559OOO+6YJk6cmDbccMN01VVX9fYuAdCPHXzwwalSqaQvfvGLc7122GGH5a/FNrXuuuuu1NLSkvbee++53vPss8/m7+npcffddxeaFihCJcuyrJBPBvqkl156Kb3yyitp4403Ti+//HLabLPN0l/+8pe01FJL9fauAdAPRTB18803p46OjryNGT58eL7+3XffTaNGjUqtra1pp512ShdffHH1Pf/yL/+SRowYkX7yk5+kxx9/PI0ePbpbwLX66qunm266Ka233nrdvmvkyJFpyJAhizF1sOj0cMEAE41fBFth5ZVXTssvv3x64403enu3AOjHNt100zR27Nh07bXXVtfF8qqrrpo22WSTbttOnz49/fd//3c69NBD8x6u2kBszuAq2qnah2CL/kjABSWzww47VIdeDB06NE2YMCH97Gc/63Hb+++/P3V2duaNJAAsis9+9rPpoosuqj6/8MIL02c+85m5trvyyivTuuuum8aPH58+/elP59sZcEWZCbigRKLBevDBB9Ppp5+eD+uIYRp77rlnOvDAA9MzzzzTbdvo1Yr1P/rRj3ptfwEojwiepk6dmp577rn8cccdd+Tr5hTDCLvWRxvV3t6ebrvttrm222abbfJhh7UP6I8G9/YOAM3zxBNPpGnTpuUNWAy9CIccckg666yz8uArxsSHGTNmpP322y99/etfzxs0AFhUK6ywQnWIYJwAjOUYtl4r2qJ77703XXfddfnzwYMHp49//ON5EBYTOtWKYYcxSgP6OwEXlEgMEVx22WXzGQjDCy+8kI499tg0bNiwfEbCEI1gXOC88847pwMOOKCX9xiAsg0rPPzww/Pl8847b67XI7CaNWtWt0kyol2Kdurcc89NbW1t1fUx3H2ttdZaTHsOxRFwQYk88MAD+dCMpZdeOr82K2aIitmiLrjggmrjFkM84qxhBGDXX399vu7SSy9NG2ywQS/vPQD9XYywmDlzZn4d8R577NHttQi0fvrTn6Yzzjgj7b777t1ei1EXl19+eY9Ty0N/J+CCkgVccc+TI444Ir311lvp6KOPTttuu223+59st912afbs2b26nwCUU9xb67HHHqsu1/rlL3+Z3nzzzXyoe21PVvjIRz6S937VBlyvv/56fvuSWssss0xaYoklCk0DNJtJM6BkAVdckxVDMDbffPN0/vnnp1NOOSW/pwkALA5x3614zCkCql133XWuYKsr4Pr973+f/vCHP1TXxbZxK5PaR9fIDOhP3PgYSuLpp59Oa665ZnrkkUfS+uuvX10f0+4edNBB6d///d97df8AAAYiPVxQogkz4oaQ66yzTrf1u+yyS3U2KAAAFi8BF5RoOOHaa6+d3+y4VgzJiGAsZiwEAGDxMqQQAACgIHq4AAAACiLgAgAAKIiACwAAoCACLgAAgIIIuAAAAAoi4AIAACiIgAsAAKAgAi4AAICCCLgAAAAKIuACAAAoiIALAACgIAIuAACAggi4AAAAUjH+P0YlTtdKsFd4AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure, axes = plt.subplots(1,2,figsize=(10,5),sharey=True)\n",
    "tukey_r2.plot_simultaneous(comparison_name=best_r2,ax=axes[0])\n",
    "r2_anova = run_anova(split_combo_df,\"R2\",\"method\")\n",
    "axes[0].set_xlabel(\"$R^2$\")\n",
    "axes[0].set_title(f\"ANOVA p={r2_anova:.2e}\")\n",
    "tukey_mae.plot_simultaneous(comparison_name=best_r2,ax=axes[1])\n",
    "axes[0].set_xlabel(\"$R^2$\")\n",
    "axes[1].set_xlabel(\"MAE\")\n",
    "mae_anova = run_anova(split_combo_df,\"MAE\",\"method\")\n",
    "axes[1].set_title(f\"ANOVA p={mae_anova:.2e}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a58fe4d8",
   "metadata": {},
   "source": [
    "Now let's use the same cross-validation function to compare machine learning algorithms.  To keep things simple we'll use scaffold splitting with three different regression models.\n",
    "- [LightGBM](https://lightgbm.readthedocs.io/en/latest/Python-Intro.html) - very fast\n",
    "- [RandomForest](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.RandomForestRegressor.html) - not so fast\n",
    "- [HistGradientBoosting](https://scikit-learn.org/stable/auto_examples/ensemble/plot_hgbt_regression.html#) - kinda fast\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14413292",
   "metadata": {},
   "source": [
    "To begin, we'll create two lists: one of regressors and another of their corresponding names.  Note that I supplied arguments to the constructors to eliminate verbose output. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "fbd43fc7-663a-4b7e-83ac-758a55ebda15",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_list = [LGBMRegressor(verbose=-1),RandomForestRegressor(),HistGradientBoostingRegressor()]\n",
    "model_names = [\"LGBM\",\"RF\",\"HGBR\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "15f07796",
   "metadata": {},
   "source": [
    "Loop over all the classifiers and compare the results of 5x5 cross-validation for each model. The model stats are stored in a dataframe and the dataframes are combined to show the comparitive results. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "411406a4-aac0-4601-8b14-ea903c08e20b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6263c6ec70934c03a49c7c79dc5f0e40",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "LGBM:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9bfb6cb14062400d862967faf1c1bc99",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "RF:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "123896f89b424d909ac2268ca2ec813e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HGBR:   0%|          | 0/5 [00:00<?, ?it/s]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "model_score_df_list = []\n",
    "for model, model_name in zip(model_list, model_names):\n",
    "    model_score_list = cross_validate(chembl_df,model,\"desc\",\"PXR_pEC50\",scaffold_clusters,model_name,metric_list)\n",
    "    model_score_df = pd.DataFrame(model_score_list,columns=metric_names)\n",
    "    model_score_df['model'] = model_name\n",
    "    model_score_df_list.append(model_score_df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e805259",
   "metadata": {},
   "source": [
    "Combine the model stats into a single dataframe. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "adaf77cb-d5c9-4eac-9b68-34f1e3e6f432",
   "metadata": {},
   "outputs": [],
   "source": [
    "model_combo_df = pd.concat(model_score_df_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c941feb5",
   "metadata": {},
   "source": [
    "Plot $R^2$ and MAE as box plots. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "5da35790-0036-4ca3-9604-c56f5dc41ee1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA90AAAHqCAYAAAAZLi26AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAATdZJREFUeJzt3Ql4VOXZ//E7CSRhX4zsgbBvyl4wolYqGqxVaWkFBIIUwYKIiqhQMKxhEUQUUxCEAgIltS91ARpLQUTKpiAvFAFlS0AgEJEEIkkgyf+6n/5n3plkAgnkzPr9XNe5knPmzJkzJOTM7zzPcz9B+fn5+QIAAAAAAEpdcOkfEgAAAAAAKEI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWKWPVgX1JXl6enD59WipVqiRBQUGePh0AQADIz8+XS5cuSZ06dSQ4mHvgt4prOQDAW6/lhG4Rc5GOjIz09GkAAALQyZMnpV69ep4+DZ/HtRwA4K3XckK3iLkrbvvHqly5sqdPBwAQADIyMkxItF2DcGu4lgMAvPVaTugWsXdD04s0F2oAgDvRFbp0cC0HAHjrtZxBZAAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAnCQkJEhUVJSEh4dLly5dZNeuXUXuu3TpUjM/qeOiz3NU8HHbMmvWLPs++noFH58xY4al7xMAAHco45ZXAQAAPiExMVFGjRolCxYsMIF77ty5EhMTI4cPH5YaNWq4fE7lypXN4zYamB2dOXPGaf0f//iHDB48WHr16uW0ffLkyTJkyBD7eqVKlUrpXQEA4DmEbgAAYDdnzhwTfAcNGmTWNXyvW7dOlixZImPGjHH5HA3ZtWrVKvKYBR/76KOPpFu3btKoUSOn7Rqyr3ccAAB8Ed3LAQCAkZOTI7t375bu3bvbtwUHB5v17du3F/m8y5cvS4MGDSQyMlIef/xxOXDgQJH7pqammhCvLd0FaXfy2267Tdq3b2+6nl+7dq0U3hUAAJ5FSzeAm5Kbmyv79u2TCxcuSPXq1aVNmzYSEhLi6dMCcAvS0tLM/+2aNWs6bdf1Q4cOuXxO8+bNTSu4/g1IT0+X2bNny913322Cd7169Qrtv2zZMtOi/Zvf/MZp+8iRI6VDhw7m78m2bdtk7Nixplu6try7kp2dbRabjIyMm3zXAABYi9ANoMS2bNkif/rTn+Ts2bP2bdoldPjw4XLfffd59NwAuFd0dLRZbDRwt2zZUt59912ZMmVKof01oPfr169QsTUdR26jAT40NFSeeeYZmT59uoSFhRU6jm6fNGlSqb8fAABKG93LAZQ4cE+YMMGMxdQKx+vXrzdfdV236+MAfFNERITpsaJdwB3penHHWpctW9Z0Dz9y5Eihx7744gtTcO3pp5++4XG0iJt2Lz9x4oTLx7UlXFvWbcvJkyeLdX4AALgboRtAsWm3U23h1latqVOnSuvWraV8+fLmq67r9vnz55v9APgebV3u2LGjbNy40b4tLy/PrDu2Zl+P/v/fv3+/1K5du9BjixcvNsdv27btDY+zd+9eM568qIrp2vqtVdMdFwAAvBHdywEUm47h1i7lr732mvkw7EjXtcvos88+a/bTli4Avke7eQ8cOFA6deoknTt3NlOGZWZm2quZx8bGSt26dU33bts0X3fddZc0adJELl68aAqgJScnF2rN1jHXH3zwgbzxxhuFXlOLtO3cudNUNNfx3rr+4osvSv/+/aVatWriL7KysiQlJUX8Xf369QsNHwCAQEboBlBsWjRNNWzY0OXjtu22/QD4nt69e8v58+clLi7O3GRr166dJCUl2YuraWh0vOn2448/minGdF8NyNqSrYXQWrVq5XTc1atXS35+vvTt29dlq7U+PnHiRFMcTf+WaOh2HOftD/TfbujQoeLvFi5cKM2aNfP0aQCA1wjK1ytggNO771WqVDFjwuieBhTt66+/Nh+EdQy3dikvSKsVa0v3m2++SUs3cANcewLv39PdLd3a4yA+Pl7GjRtnpnRzF1q6AQSKjGJee2jpBlBsWlFYiymtXLnSjOF2bO3ScZ+6Xcdx6n4AAGcaRD3RAqyBm5ZnAPAcCqkBKDataqzTgul4y/Hjx5uW7Z9++sl81XXdPmzYMObrBgAAAP4/WroBlIjOw61z42oVc+1KbqMt3LqdeboBAACA/0PoBlBiGqy7du1qqpRr0bTq1aubLuW0cAMAAADOCN0AbooGbIqlAQAAANfHmG4AAAAAACxC6AYAAAAAwCJ0LwcAAAAAH5OVlSUpKSniz+rXr2+mW/R1hG4AAAAA8DEauIcOHSr+bOHChdKsWTPxdYRuAAAAAPDBVmANpe6SnJws8fHxMm7cOGnQoIHb3qM/IHT7iEDoPuJPXUgAAAAAK+lnZk+0Amvg9ofWZ3cidPuIQOg+4k9dSAAAAABAEbp9RCB0H/GnLiQAAAAAoAjdPoLuIwAAAADge7xynu6EhASJiooyQbNLly6ya9eu6+5/8eJFefbZZ6V27doSFhZmQuL69evddr4AAAAAAPhES3diYqKMGjVKFixYYAL33LlzJSYmRg4fPiw1atQotH9OTo48+OCD5rG//e1vUrduXdM1umrVqh45fwAAAAAAvDZ0z5kzR4YMGSKDBg0y6xq+161bJ0uWLJExY8YU2l+3X7hwQbZt2yZly5Y127SVHAAAAAAAT/Oq0K2t1rt375axY8fatwUHB0v37t1l+/btLp/z8ccfS3R0tOle/tFHH8ntt98uTz75pLz66qsSEhLi8jnZ2dlmscnIyLDg3QCA92IaQgAAgAAM3WlpaZKbmys1a9Z02q7rhw4dcvmcY8eOyaZNm6Rfv35mHPeRI0dk+PDhcvXqVZkwYYLL50yfPl0mTZpkyXsAAF/ANIQAAAABGLpvRl5enhnPrR+stGW7Y8eO8v3338usWbOKDN3akq7jxh1buiMjI9141gDgWUxDCAAAEIChOyIiwgTn1NRUp+26XqtWLZfP0YrlOpbbsSt5y5Yt5ezZs6a7emhoaKHnaIVzXQAgUDENIQAAQABOGaYBWVuqN27c6NSSres6btuVrl27mi7lup/Nt99+a8K4q8ANAAAAAEBAhm6l3b4XLVoky5Ytk4MHD8qwYcMkMzPTXs08NjbWqdCaPq7Vy59//nkTtrXS+bRp00xhNQAAAAAAPMmruper3r17y/nz5yUuLs50EW/Xrp0kJSXZi6tp8R+taG6jY7E//fRTefHFF6VNmzZmnm4N4Fq9HAAAAAAAT/K60K1GjBhhFlc2b95caJt2Pd+xY4cbzgwAAAAAAB/uXg4AAAAAgL8gdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAcJKQkCBRUVESHh4uXbp0kV27dhW579KlSyUoKMhp0ec5euqppwrt06NHD6d9Lly4IP369ZPKlStL1apVZfDgwXL58mXL3iMAAO5Sxm2vBAAAvF5iYqKMGjVKFixYYAL33LlzJSYmRg4fPiw1atRw+RwNyvq4jYbqgjRk//nPf7avh4WFOT2ugfvMmTOyYcMGuXr1qgwaNEiGDh0qq1atKtX3BwCAuxG6AQCA3Zw5c2TIkCEm9CoN3+vWrZMlS5bImDFjXD5HQ3atWrWue1wN2UXtc/DgQUlKSpIvv/xSOnXqZLbNmzdPfvnLX8rs2bOlTp06t/y+AADwFLqXAwAAIycnR3bv3i3du3e3bwsODjbr27dvL/J52g28QYMGEhkZKY8//rgcOHCg0D6bN282LeXNmzeXYcOGyQ8//GB/TI+tXcptgVvpa+pr79y5s1TfIwAA7kboBgAARlpamuTm5krNmjWdtuv62bNnXT5HQ7S2gn/00UeyYsUKycvLk7vvvltOnTrl1LV8+fLlsnHjRpk5c6Z8/vnn8vDDD5vXUnrsgl3Xy5QpI9WrVy/ydbOzsyUjI8NpAQDAG9G9HAAA3LTo6Giz2Gjgbtmypbz77rsyZcoUs61Pnz72x++8805p06aNNG7c2LR+P/DAAzf1utOnT5dJkyaVwjsAAMBatHQDAAAjIiJCQkJCJDU11Wm7rt9ozLZN2bJlpX379nLkyJEi92nUqJF5Lds+euxz58457XPt2jVT0byo1x07dqykp6fbl5MnTxbr/AAAcDdCNwAAMEJDQ6Vjx46mG7iNdhfXdcfW7OvRLuP79++X2rVrF7mPdj3XMd22ffTYFy9eNOPJbTZt2mReWyuoF1WYTaumOy4AAHgjQjcAALDT6cIWLVoky5YtM1XFtehZZmamvZp5bGysaWW2mTx5svzzn/+UY8eOyZ49e6R///6SnJwsTz/9tL3I2ssvvyw7duyQEydOmACvxdaaNGlipiJT2h1dx31r1XSdE/zf//63jBgxwnRLp3I5AMDXMaYbAADY9e7dW86fPy9xcXGmiFm7du3MdF624mopKSmmqrjNjz/+aMKy7lutWjXTUr5t2zZp1aqVeVy7q+/bt8+EeG3N1hD90EMPmfHejnN1r1y50gRtHeOtx+/Vq5e8/fbbHvgXAACgdBG6AQCAEw2/uriixc8cvfnmm2YpSrly5eTTTz+94WtqpfJVq1bdxNkCAODd6F4OAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUY0w34maysLFPoyJ/Vr19fwsPDPX0aAAAAwA0RugE/o4F76NCh4s8WLlwozZo18/RpAAAAADdE6Ab8sBVYQ6k76Fy88fHxMm7cOGnQoIG48z0CAAAAvoDQDfgZ7Xbt7lZgDdy0PAMAAACFUUgNAAAAAACLELoBAAAAALAIoRsAAAAAAIsQugEAAAAAsAihGwAAAAAAixC6AQAAAACwCKEbAAAAAACLELoBAAAAALAIoRsAAAAAAIsQugEAAAAAsEgZqw4MAADgC1JTUyU9PV38TXJystNXf1OlShWpWbOmp08DAG6I0A0AAAI6cPcfECtXc7LFX8XHx4s/KhsaJiveX07wBuD1CN0AACBgaQu3Bu4rjX4ueeFVPH06KKbgrHSRY5+bnx+hG4C3I3QDAICAp4E7r0KEp08DAOCHCN0AAACAF8nKypKUlBTxZ/Xr15fw8HBPnwbgFoRuAAAAwIto4B46dKj4s4ULF0qzZs08fRqAWxC6AQAAAC9rBdZQ6g5a3V6L7Y0bN04aNGgg7nyPQKAgdAMAAABeRLtdu7sVWAM3Lc+ANYItOi4AAAAAAAGP0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAARS6E5ISJCoqCgzR2GXLl1k165dRe67dOlSCQoKclr0eQAAAAAAeFoZ8TKJiYkyatQoWbBggQncc+fOlZiYGDl8+LDUqFHD5XMqV65sHrfR4O0uqampkp6eLv4mOTnZ6au/qVKlitSsWdPTpwEAAADAz3ld6J4zZ44MGTJEBg0aZNY1fK9bt06WLFkiY8aMcfkcDdm1atVy85n+N3D3HxArV3OyxV/Fx8eLPyobGiYr3l9O8IbX4Uaeb+JGHgAA8InQnZOTI7t375axY8fatwUHB0v37t1l+/btRT7v8uXL0qBBA8nLy5MOHTrItGnTpHXr1kXun52dbRabjIyMmzpf/WCsgftKo59LXniVmzoG3C84K13k2Ofm58eHZHgTbuT5Lm7kAQAAnwjdaWlpkpubW+hDi64fOnTI5XOaN29uWsHbtGljQtTs2bPl7rvvlgMHDki9evVcPmf69OkyadKkUjtvDdx5FSJK7XgAAhM38nwTN/IAAP7eY03Ra81PQvfNiI6ONouNBu6WLVvKu+++K1OmTHH5HG1J13Hjji3dkZGRbjlfALgRbuQBAOB7AqHHmqLXmo+H7oiICAkJCTG/sI50vbhjtsuWLSvt27eXI0eOFLlPWFiYWQAAAACgNNBjzXcFW9xrzatCd2hoqHTs2FE2btwoPXv2NNt0nLaujxgxoljH0O7p+/fvl1/+8pcWny0AAP5Jp+6cNWuWnD17Vtq2bSvz5s2Tzp07Fzl1p634qY3e2M7KyjLfX716VcaPHy/r16+XY8eOme57WqtlxowZUqdOHftzdKrQgl0WdThYUUVUAcBb0WMNXh26lXb7HjhwoHTq1Mlc4HXKsMzMTPsFPTY2VurWrWsuxGry5Mly1113SZMmTeTixYvmQ4JetJ9++mkPvxMAAHxPaU/d+dNPP8mePXvktddeMwH+xx9/lOeff14ee+wx+eqrr5yOo9d0ncHEplKlSpa8RwAAAjp09+7dW86fPy9xcXHmDnu7du0kKSnJ3syfkpJiKprb6MVbL9C6b7Vq1UxL+bZt26RVq1YefBcAAPim0p66U1u2N2zY4LTtnXfeMTfW9Zpev359p5DtiSlAAQCw0v+lVy+iXcm1tVqn9dq5c6e5026zefNm05XN5s0337Tvq8FbPxjomG4AAHBzU3dq9++bmbpTi5I+/vjjZgaR69ExcxrUq1at6rRdu5zfdttt5jquPdeuXbtWCu8KAADP8rqWbsAf+ev0Ef4+dYTV00cA3sYdU3fqWO9XX31V+vbta7ql24wcOVI6dOgg1atXNz3WdKaRM2fOmJZ3V/Rmuy6OM5EAAOCNCN2AxQJh+gh/nTrC6ukjAH9Qkqk7tajaE088Ifn5+TJ//nynxxyn8tQAr8VVn3nmGVPDxdWMI7p90qRJlrwnAABKE6EbsBjTR/guq6ePALyNlVN32gK39ozZtGmTUyu3Kzq0TLuXnzhxwrSmF6Qt4Y5BXVu6tXs7AADehtANuAnTRwDwdlZN3WkL3N9995189tlnZtz2jezdu9eMJy+qYrq2frtqAQcAwNsQugEAgGVTd2rg/u1vf2umDVu7dq0J5Vr4VOn4bQ36WqRNC6d269bNVDDX9RdffFH69+9vZiYBAMCXEboBAIBlU3d+//338vHHH5vv9ViOtNX7/vvvNy3Wq1evlokTJ5riaA0bNjSh27H7OAAAvorQDQAAnGhX8qK6k+vUnY506k5dihIVFWUKp12PVi3fsWPHTZ4tAADezSvn6QYAAAAAwB8QugEAAAAAsAihGwAAAAAAixC6AQAAAACwCKEbAAAAAACLELoBAAAAALAIoRsAAAAAAIsQugEAAAAAsAihGwAAAAAAi5Sx6sAAnAVfuejpU0AJ8TMDAADArSJ0A25S7vgWT58CAAAAADcjdANucqXhfZJXrqqnTwMlbOnmZgkAAABuBaEbcBMN3HkVIjx9GgAAAADciEJqAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhHm6AcDLBF+56OlTQAnw8wIAANdD6AYAL1Pu+BZPnwIAAABKCaEbALzMlYb3SV65qp4+DZSgpZsbJQAAoCiEbgDwMhq48ypEePo0AAAAUAoopAYAAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFiE0A0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdAAAAAABYhNANAAAAAIBFCN0AAAAAAFikjFUHBgAAAPxFamqqpKeni79JTk52+uqPqlSpIjVr1vT0aSCAEboBAACAGwTu/gNi5WpOtvir+Ph48VdlQ8NkxfvLCd7wGEI3AAAAcB3awq2B+0qjn0teeBVPnw5KIDgrXeTY5+ZnSOiGpxC6AQCAk4SEBJk1a5acPXtW2rZtK/PmzZPOnTu73Hfp0qUyaNAgp21hYWGSlZVlX8/Pz5cJEybIokWL5OLFi9K1a1eZP3++NG3a1L7PhQsX5LnnnpNPPvlEgoODpVevXvLWW29JxYoVLXynQMlo4M6rEOHp0wDgYyikBgAA7BITE2XUqFEmJO/Zs8eE7piYGDl37lyRz6lcubKcOXPGvhQcG/r666/L22+/LQsWLJCdO3dKhQoVzDEdg3m/fv3kwIEDsmHDBlm7dq1s2bJFhg4daul7BQDAHQjdAADAbs6cOTJkyBDTet2qVSsTlMuXLy9Lliwp8jlBQUFSq1Yt++LYhVNbuefOnSvjx4+Xxx9/XNq0aSPLly+X06dPy4cffmj2OXjwoCQlJcl7770nXbp0kXvuuce0rq9evdrsBwCALyN0AwAAIycnR3bv3i3du3e3b9Ou3rq+ffv2Ip93+fJladCggURGRppgrS3WNsePHzfd1B2PqZWENVzbjqlfq1atKp06dbLvo/vra2vLuCvZ2dmSkZHhtAAA4I0I3QAAwEhLS5Pc3NxCxYZ0XYOzK82bNzet4B999JGsWLFC8vLy5O6775ZTp06Zx23Pu94x9WuNGjWcHi9TpoxUr169yNedPn26Ce+2RQM/AADeiNANAABuWnR0tMTGxkq7du3k5z//uaxZs0Zuv/12effddy193bFjx5pqxLbl5MmTlr4eAAA3i9ANAACMiIgICQkJMXMSO9J1HatdHGXLlpX27dvLkSNHzLrtedc7pn4tWKjt2rVrpqJ5Ua+rFdK1gJvjAgCANyJ0AwAAIzQ0VDp27CgbN260b9Pu4rquLdrFod3T9+/fL7Vr1zbrDRs2NMHZ8Zg6/lrHatuOqV91KjEdT26zadMm89o69hsAAF/GPN0AAMBOpwsbOHCgKWqmc3Nr5fHMzEz7XNzalbxu3bpmTLWaPHmy3HXXXdKkSRMTnHV+b50y7Omnn7ZXNn/hhRdk6tSpZl5uDeGvvfaa1KlTR3r27Gn2admypfTo0cNUTddq6VevXpURI0ZInz59zH4AAPgyQjcAALDr3bu3nD9/XuLi4kwRMx2rrdN52QqhpaSkmKriNj/++KMJy7pvtWrVTEv5tm3bzHRjNq+88ooJ7jrvtgZznRJMjxkeHm7fZ+XKlSZoP/DAA+b4vXr1MnN7AwDg67yye3lCQoJERUWZi7F2K9u1a1exnqfzeeodddudcwAAUHIafrW1Wqfl0m7gjl28N2/eLEuXLrWvv/nmm/Z9NXivW7fOjOl2pNdmbRHXx7OysuRf//qXNGvWzGkfrVS+atUquXTpkimMphXRK1as6IZ3CwBAgLV0JyYmmq5t2r1ML/LarS0mJkYOHz5caDoRRydOnJDRo0fLvffe69bzBQAAvi/4ykVPnwJKgJ8XAF/idaF7zpw5ppuabeyYhm+9a653vMeMGVNk0ZZ+/frJpEmT5IsvvjBd1wAAAIqr3PEtnj4FAICf8qrQnZOTYyqX6tybNjquq3v37rJ9+/Yin6dd1rQVfPDgwSZ0AwAAlMSVhvdJXrmqnj4NlKClmxslAHyFV4XutLQ002ptK9Zio+uHDh1y+ZytW7fK4sWLZe/evcV+HR13povj1CUAACBwaeDOqxDh6dMAAPghrwrdJaXFVgYMGCCLFi2SiIjiXyh1mhPtig64U3BWuqdPASXEzwwAAAB+Fbo1OIeEhEhqaqrTdl2vVatWof2PHj1qCqg9+uij9m15eXnma5kyZUzxtcaNGxd6nnZf12Jtji3dkZGRpfxugP+qUqWKlA0NEzn2uadPBTdBf3b6MwQAAAB8PnSHhoaa+T03btxon/ZLQ7Su6/QlBbVo0UL279/vtG38+PGmBfytt94qMkiHhYWZBXAHHR6x4v3lZgocf6PTBMXHx8u4ceOkQYMG4o80cBcc8gIAAAD4ZOhW2gI9cOBA6dSpk3Tu3NlMGZaZmWmvZh4bGyt169Y1XcR1Hu877rjD6flVq/63CErB7YAnaWjz5+CmgbvgnLsAAAAAvDB09+7dW86fPy9xcXFy9uxZadeunSQlJdkDS0pKiqloDgAAAACAt/O60K20K7mr7uRq8+bN133u0qVLLTorAAAAAABKhiZjAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxSxqoDAwAAAECgCb5y0dOnAC/7mRG6AQAAAKCUlDu+xdOnAC9D6AYAAACAUnKl4X2SV66qp08DJWzptvJmCaEbAAAAAEqJBu68ChGePg14EQqpAQAAAABgEUI3AAA+ateuXZKbm1vk49nZ2fLXv/7VrecEAACcEboBAPBR0dHR8sMPP9jXK1euLMeOHbOvX7x4Ufr27euhswMAAIrQDQCAj8rPz7/uelHbAACA+xC6AQDwY0FBQZ4+BQAAAhqhGwAAAAAAb5ky7MqVK3LhwgWpW7eu0/YDBw5I69atS/PcAADADXzzzTdy9uxZe1fyQ4cOyeXLl816Wlqah88O8L+5fOFb+JnB50L33/72N3nhhRckIiJC8vLyZNGiRdKlSxfz2IABA2TPnj1WnScAAHDhgQcecBq3/atf/crerVy3070cKD3ljm/x9CkA8PfQPXXqVNm9e7fUrFnTfB04cKD88Y9/lCeffJJCLQAAuNnx48c9fQpAQLnS8D7JK1fV06eBErZ0c7MEPhW6r169agK36tixo2zZskV+/etfy5EjR7iTDgCAmzVo0OCG+/znP/9xy7kAgUADd16FCE+fBgB/LqRWo0YN2bdvn329evXqsmHDBjl48KDTdgAA4DmXLl2ShQsXSufOnaVt27Ylfn5CQoJERUVJeHi4GUa2a9euYj1v9erV5iZ8z549nbbrNlfLrFmz7Pvo6xV8fMaMGSU+dwAAfDp0v//++yZ4OwoNDZW//OUv8vnnn5f2uQEAgBLQHmg69Kt27doye/Zs+cUvfiE7duwo0TESExNl1KhRMmHCBFOrRUN7TEyMnDt37rrPO3HihIwePVruvffeQo+dOXPGaVmyZIkJ1b169XLab/LkyU77PffccyU6dwAAfL57eb169Yp8rGvXrqVxPgAQ8IKz0j19CvChn5dWLl+6dKksXrxYMjIy5IknnpDs7Gz58MMPpVWrViU+3pw5c2TIkCEyaNAgs75gwQJZt26dCcpjxoxx+Zzc3Fzp16+fTJo0Sb744gu5eNG5WnCtWrWc1j/66CPp1q2bNGrUyGl7pUqVCu0LAEDATRnmKDk5WQ4fPixt2rRxeZE8ffq01KlT51ZeAgACRpUqVaRsaJjIMXoO+Rr9uenPz90effRR07r9yCOPyNy5c6VHjx4SEhJigvLNyMnJMYVSx44da98WHBws3bt3l+3btxf5PG2h1p5wgwcPNqH7elJTU02IX7ZsWaHHtDv5lClTpH79+qZI64svvihlyrj+qKI3FnSx0RsOAAD4VejWLuWxsbHm7raO+Xr33XfNtGEpKSmyatUq+fvf/24u3NeuXSvdMwYAP6WFKle8v1zS0/2vpVtv0sbHx8u4ceOKVfzL12jgthUadad//OMfMnLkSBk2bJg0bdr0lo+n83rrdb3ge9F1nf/bla1bt5pW9r179xbrNTRsa4v2b37zG6ft+j46dOhg6sVs27bNBH/tYq4t765Mnz7dtKwDAOC3oVvvROtYK72rrdOG6QX/22+/NXepGzdubOYNLaobGgDANQ03nghv7qKBu1mzZp4+Db9hC7w6o0jLli3Nze8+ffq4tWCbvuaiRYskIqJ4FZ21m7p2Rdcb9o50HLmN9qDTmjHPPPOMCddhYWGFjqOh3PE52tIdGRl5S+8HAACvCt1Hjx6V559/3nyA0iqn2hXs3//+t6lirhf+QJv/D76DnxcAf3HXXXeZRbuWawE0DbQaRPPy8szsIhpCtVW5uDQ4a/d07QLuSNddDSPTzwJaQE27udvoayvtFq5D0PRGvI12Pddteq43olXTtbecHr958+aFHtcg7iqMAwDgN6Fb5+wuV66cvcCa3rHWSqmBFrhVueNbPH0KAIAAVqFCBfn9739vFg212vqtPc+0x9mDDz4oH3/8cbGOo63L2mq+ceNG+7RfGqJ1fcSIEYX2b9Gihezfv99p2/jx400L+FtvvVWo5dnWKl+cacy0u7qOJy84awoAAAFVSE3HbmvRFr3o6p3xatWqSSC60vA+yStX1dOngRK0dHOjBIC/0lbh119/3XTLXrt2rWn9LgltKddpxzp16mTm+dZW9MzMTHs1c63nUrduXXN8veF+xx13OD2/atX/Xg8Lbtfu3x988IG88cYbhV5Ti7Tt3LnTVDTXlnld1yJq/fv3D9jPFgAA/3HToVvn4dQ5PF966SVzQczKyjJ3te+++25zodUxe0VVHPU3GrjzKhRvLBsAAKVFW7Zv5LbbbivRMXv37i3nz5+XuLg4Mx1Zu3btJCkpyV5rQAumagt0Sa1evVry8/Olb9++hR7TbuL6+MSJE01F8oYNG5rQ7ThmGwAAX3XTqfjzz/87pc13331nqpTv2bPHLMuXLzfzc2oXNQ3eOsYbAACUPp2fW2urtG/f3gRaV4KCgkp8XO1K7qo7udq8efMNz8mVoUOHmsUVrVq+Y8eOEp8nAAC+4JabonWKEl0cq6UeP35cvvrqK/n6669v9fAAAKAIOnOITuGp113t/q3dsXXKLQAA4D1K3j+sGLRb2O9+9zuZNm2aFYcHAAAiZvYQncv6lVdekU8++cQULnviiSfk008/LbLlGwAA+EHoBgAA7qHjoXWctE4R9s0330jr1q1l+PDhEhUVJZcvX/b06QEAEPAI3QAA+AktcKZjuLWVOzc319OnAwAACN0AAPg2rfat47p1Pm4tYKrzZr/zzjumynjFihU9fXoAAAS8wJjTCwAAP6TdyHWqLR3LrdOHafiOiGAKSwAAvAmhGwAAH7VgwQKpX7++NGrUyEzlaZvOs6A1a9a4/dwAAMB/EboBAPBRsbGxNzUPNwAAcB9CNwAAPmrp0qWePgUAAHADFFIDAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIlQvBwAAAS84K93Tp4AS4OcFwJcQugEAQMCqUqWKlA0NEzn2uadPBSWkPzf9+QGAtyN0AwCAgFWzZk1Z8f5ySU/3v5bT5ORkiY+Pl3HjxkmDBg3E32jg1p8fAHg7rwzdCQkJMmvWLDl79qy0bdtW5s2bJ507d3a575o1a2TatGly5MgRuXr1qjRt2lReeuklGTBggNvPG/AGWVlZkpKS4rYPdI5f3aV+/foSHh7u1tcE4L80uPlzeNPA3axZM0+fBgAELK8L3YmJiTJq1ChZsGCBdOnSRebOnSsxMTFy+PBhqVGjRqH9q1evbu7gtmjRQkJDQ2Xt2rUyaNAgs68+Dwg0GriHDh3q1tfUlhR3WrhwIR8gAQAA4BO8LnTPmTNHhgwZYoKz0vC9bt06WbJkiYwZM6bQ/vfff7/T+vPPPy/Lli2TrVu3EroRkLQVWEOpv79HAAAAwBd4VejOycmR3bt3y9ixY+3bgoODpXv37rJ9+/YbPj8/P182bdpkWsVnzpxZ5H7Z2dlmscnIyCiFswe8g3a7phUYAAAA8A5eNU93Wlqa5ObmFhpXpes6vrsoWvykYsWKpnv5I488YsaAP/jgg0XuP336dFN8w7ZERkaW6vsAAAAAAMDrQvfNqlSpkuzdu1e+/PJLM7ZUx4Rv3ry5yP21JV2Dum05efKkW88XAAAAABAYvKp7eUREhISEhEhqaqrTdl2vVatWkc/TLuhNmjQx37dr104OHjxoWrMLjve2CQsLMwsAAAAAAAHT0q3dwzt27CgbN260b8vLyzPr0dHRxT6OPsdxzDYAAAAAABLoLd1Ku4YPHDhQOnXqZObm1inDMjMz7dXMY2NjpW7duqYlW+lX3bdx48YmaK9fv17ef/99mT9/voffCQAAAAAg0Hld6O7du7ecP39e4uLiTPE07S6elJRkL66mcxBrd3IbDeTDhw+XU6dOSbly5cx83StWrDDHAQAAAADAk7wudKsRI0aYxZWCBdKmTp1qFgAAAAAAvI1Xhm4AAAAA8EXBWemePgV42c+M0A0AAAAAt6hKlSpSNjRM5Njnnj4V3AT92enP0AqEbgAAAAC4RVqDasX7yyU93T9bupOTkyU+Pl7GjRsnDRo0EH9TpUoVex2x0kboBgAAAIqBbsO+x90/Mw1tVgU3b6GBu1mzZp4+DZ9C6AYAAACug27Dvs3KbsNAcRC6AQAAgADtNuzvXYat7jYMFAehGwAAAAjwbsN0GQasQ+guBYzv8S38vAAAAAC4C6H7FjC+x3cxtgcAAACAOxC6bwHje3wXY3sAAAAAuAOh+xYxvgcAAAAAUJTgIh8BAAABKSEhQaKioiQ8PFy6dOkiu3btKtbzVq9eLUFBQdKzZ0+n7U899ZTZ7rj06NHDaZ8LFy5Iv379pHLlylK1alUZPHiwXL58uVTfFwAAnkDoBgAAdomJiTJq1CiZMGGC7NmzR9q2bSsxMTFy7ty56z7vxIkTMnr0aLn33ntdPq4h+8yZM/blL3/5i9PjGrgPHDggGzZskLVr18qWLVtk6NChpfreAADwBEI3AACwmzNnjgwZMkQGDRokrVq1kgULFkj58uVlyZIlRT4nNzfXhOZJkyZJo0aNXO4TFhYmtWrVsi/VqlWzP3bw4EFJSkqS9957z7Ss33PPPTJv3jzTcn769GlL3icAAO5C6AYAAEZOTo7s3r1bunfvbt8WHBxs1rdv317k8yZPniw1atQwXcKLsnnzZrNP8+bNZdiwYfLDDz/YH9Nja5fyTp062bfpa+pr79y50+XxsrOzJSMjw2kBAMAbEboBAICRlpZmWq0LFgjV9bNnz7p8ztatW2Xx4sWyaNGiIo+rXcuXL18uGzdulJkzZ8rnn38uDz/8sHktpcfWQO6oTJkyUr169SJfd/r06WYmCtsSGRl5E+8YAADrUb0cAADclEuXLsmAAQNM4I6IiChyvz59+ti/v/POO6VNmzbSuHFj0/r9wAMP3NRrjx071ow9t9GWboI3AMAbEboBAIChwTkkJERSU1Odtuu6jsMu6OjRo6aA2qOPPmrflpeXZ2+pPnz4sAnXBem4b32tI0eOmNCtxy5YqO3atWumormr17WNEdcFAABvR/dyAABghIaGSseOHU03cMcQrevR0dGF9m/RooXs379f9u7da18ee+wx6datm/m+qJbnU6dOmTHdtWvXNut67IsXL5rx5DabNm0yr62F1QAA8GW0dAMAADvtsj1w4EBT1Kxz584yd+5cyczMNNXMVWxsrNStW9eMqdZ5vO+44w6n52tBNGXbrnNta1XzXr16mVZrbR1/5ZVXpEmTJmYqMtWyZUsz7lurpmu19KtXr8qIESNMt/Q6deq4/d8AAIDSROgGAAB2vXv3lvPnz0tcXJwpYtauXTsznZetuFpKSoqpKl5c2l193759smzZMtOarSH6oYcekilTpjh1D1+5cqUJ2trdXI+vIf3tt9+25D0CAOBOhG4AAOBEw68urmjxs+tZunSp03q5cuXk008/veFraqXyVatWlfBMAQDwfozpBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAAAIpdCckJEhUVJSEh4dLly5dZNeuXUXuu2jRIrn33nulWrVqZunevft19wcAAAAAIGBDd2JioowaNUomTJgge/bskbZt20pMTIycO3fO5f6bN2+Wvn37ymeffSbbt2+XyMhIeeihh+T77793+7kDAAAAAODVoXvOnDkyZMgQGTRokLRq1UoWLFgg5cuXlyVLlrjcf+XKlTJ8+HBp166dtGjRQt577z3Jy8uTjRs3uv3cgUCSm5srX3/9tfm/pl91HQAAAICzMuJFcnJyZPfu3TJ27Fj7tuDgYNNlXFuxi+Onn36Sq1evSvXq1S08UyCwbdmyRf70pz/J2bNn7dtq1aplboDdd999Hj03AAAAwJt4VUt3WlqaaS2rWbOm03Zdd/xwfz2vvvqq1KlTxwT1omRnZ0tGRobTAqD4gVuHfzRq1MjUX1i/fr35quu6XR8H4NtKUlvF0erVqyUoKEh69uxp36Y3wvXafOedd0qFChXMNTo2NlZOnz7t9Fx9PX2u4zJjxoxSf28AAAR06L5VenHWC/7f//5380GhKNOnT5cqVarYFx0HDuDG9KaYtnBHR0fL1KlTpXXr1mb4h37Vdd0+f/58upoDPqyktVVsTpw4IaNHjzbFTQv2QNPjvPbaa+brmjVr5PDhw/LYY48VOsbkyZPlzJkz9uW5554r9fcHAEBAh+6IiAgJCQmR1NRUp+26rl1Xr2f27NkmdP/zn/+UNm3aXHdf7b6enp5uX06ePFkq5w/4u3379pleJ/369TNDPxzpum7XD8q6HwDfVNLaKkpvtOn//0mTJpleL4705vaGDRvkiSeekObNm8tdd90l77zzjhlOlpKS4rRvpUqVzPXetmjLOAAAvs6rQndoaKh07NjRqQiarSiatqAV5fXXX5cpU6ZIUlKSdOrU6YavExYWJpUrV3ZaANzYhQsXzNeGDRu6fNy23bYfAN9iq63iOESrOLVVtIW6Ro0aMnjw4GK9jt7w1u7jVatWddquN89vu+02ad++vcyaNUuuXbt2C+8GAADv4FWF1JR2aRs4cKAJz507d5a5c+dKZmamueOudBxY3bp1TRdxNXPmTImLi5NVq1aZ8WC2sd8VK1Y0C4DSYytQePz4cdOlvCDd7rgfAN9yvdoqhw4dcvmcrVu3yuLFi2Xv3r3Feo2srCwzxlun+3S86T1y5Ejp0KGD+fuxbds20ytNe85oy3tR9Vl0saE+CwDAW3ld6O7du7ecP3/eBGkN0DoVmLZg2z4AaFc0x26tOn5U78z/9re/dTqOjkWbOHGi288f8Gc6dEO7fOpUfTqG2/H/ovZK0e21a9e+4RAPAP7h0qVLMmDAAFm0aJEZInYjWlRNu5nn5+eb63fBm+42+jdEe78988wz5ia79lArSLdrd3YAALyd14VuNWLECLO4snnz5kKFWwC4h9Zc0GnB9KbW+PHjzRhO7VKuLdwauLX7qX4I1v0A+J6S1lY5evSouQ4/+uijTjfgVJkyZUzBtMaNGzsF7uTkZNm0adMNh3Zp1XTtXq7H17HgBWlLuGNQ15ZuCqMCALyRV4ZuAN5L5+HWYK1VzJ999ln7dm3h1u3M0w34LsfaKrZpv2y1VVzdDG/RooXs37/faZvekNMW8Lfeessegm2B+7vvvpPPPvvMjNu+Ee2urr1pdKy4K9r67aoFHAAAb0PoBlBiGqy7du1qqpRr0TQdg6ndQWnhBnxfSWqr6PScd9xxh9PzbcXRbNs1cOsQMJ0ubO3atWbMuK3+iv7t0KCvvWR27twp3bp1MxXMdf3FF1+U/v37S7Vq1dz+bwAAQGkidAO4KRqwtcIwAP9S0toqN/L999/Lxx9/bL7XYznSVu/777/ftFivXr3a1GLR4mg6bEVDt2P3cQAAfBWhGwAA3HRtlYKWLl3qtK4zi2jhtOvRquU7duy4iTMFAMD7edU83QAAAAAA+BNCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYpY9WBAQDeKysrS1JSUtz2esnJyU5f3aV+/foSHh7u1tcEAABwROgGgACkgXvo0KFuf934+Hi3vt7ChQulWbNmbn1NAAAAR4RuAAhA2gKsgTQQ3icAAIAnEboBIABpl2tagAEAAKxHITUAAAAAACxCSzcAAAAA+JhAKIpa308KohK6AQAAAMDHBEJR1IV+UhCV0A0AAAAAPiYQiqLW95OCqIRuAAAAAPAxFEX1HRRSAwAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsQugGAAAAAMAihG4AAAAAACxC6AYAAAAAwCKEbgAAAAAALELoBgAAAADAIoRuAAAAAAAsUsaqAwMAAAAouaysLElJSXHLayUnJzt9dZf69etLeHi4W18T8BRCNwAAAOBFNHAPHTrUra8ZHx/v1tdbuHChNGvWzK2vCXgKoRsAAADwItoKrKHU398jECgI3QAAAIAX0W7XtAID/oNCagAAwElCQoJERUWZD/5dunSRXbt2Fet5q1evlqCgIOnZs6fT9vz8fImLi5PatWtLuXLlpHv37vLdd9857XPhwgXp16+fVK5cWapWrSqDBw+Wy5cvl+r7AgDAEwjdAADALjExUUaNGiUTJkyQPXv2SNu2bSUmJkbOnTt33eedOHFCRo8eLffee2+hx15//XV5++23ZcGCBbJz506pUKGCOaYWi7LRwH3gwAHZsGGDrF27VrZs2eL2Ma0AAFiB0A0AAOzmzJkjQ4YMkUGDBkmrVq1MUC5fvrwsWbKkyOfk5uaa0Dxp0iRp1KhRoVbuuXPnyvjx4+Xxxx+XNm3ayPLly+X06dPy4Ycfmn0OHjwoSUlJ8t5775mW9XvuuUfmzZtnWs51PwAAfBljun2EO6eOUEwfAQCBJycnR3bv3i1jx461bwsODjbdwbdv317k8yZPniw1atQwXcK/+OILp8eOHz8uZ8+eNcewqVKlignXesw+ffqYr9qlvFOnTvZ9dH99bW0Z//Wvf13q7xUAAHchdPsIT0wdoZg+AgACR1pammm1rlmzptN2XT906JDL52zdulUWL14se/fudfm4Bm7bMQoe0/aYftXQ7qhMmTJSvXp1+z4FZWdnm8UmIyOjWO8RAAB3I3T7iECYOkIxfQQA+I5Lly7JgAEDZNGiRRIREeHW154+fbrpzg4AgLcjdPsIpo4AAFhNg3NISIikpqY6bdf1WrVqFdr/6NGjpoDao48+at+Wl5dnb6k+fPiw/Xl6DK1e7njMdu3ame91n4KF2q5du2Yqmrt6XaVd4LXgm2NLd2Rk5E2+cwAArEMhNQAAYISGhkrHjh1l48aNTiFa16Ojowvt36JFC9m/f7/pWm5bHnvsMenWrZv5XkNww4YNTXB2PKYGZB2rbTumfr148aIZT26zadMm89o69tuVsLAwM72Y4wIAgDeipRsAANhp6/HAgQNNUbPOnTubyuOZmZmmmrmKjY2VunXrmu7d2gvrjjvucHq+FkRTjttfeOEFmTp1qjRt2tSE8Ndee03q1Kljn8+7ZcuW0qNHD1M1XaulX716VUaMGGGKrOl+AAD4MkI3AACw6927t5w/f17i4uJMETPtAq7TedkKoWlhT60qXhKvvPKKCe5aEFRbtHVKMD2m42wVK1euNEH7gQceMMfv1auXmdsbAABfF5SvE2gGOO3mptOXpKen0z0NAOAWXHtKF/+ehX377bfmRgczgwCAZ689jOkGAAAAAMAihG4AAAAAACxC6AYAAAAAIJBCd0JCgkRFRZkCKzpVyK5du4rc98CBA6bYiu4fFBRkqqwCAAAAAOANvC50JyYmmulKJkyYIHv27JG2bdtKTEyMnDt3zuX+P/30kzRq1EhmzJhh5gEFAAAAAMBbeF3onjNnjpmnU+cDbdWqlZmvs3z58rJkyRKX+//sZz+TWbNmmbk8w8LC3H6+AAAAAAD4ROjOycmR3bt3S/fu3e3bdK5OXd++fXupvU52drYp7+64AAAAAADg16E7LS1NcnNzpWbNmk7bdf3s2bOl9jrTp08386nZlsjIyFI7NgAAAAAAXhm63WXs2LFmAnPbcvLkSU+fEgAAAADAD5URLxIRESEhISGSmprqtF3XS7NImo79Zvw3AAAAACCgWrpDQ0OlY8eOsnHjRvu2vLw8sx4dHe3RcwMAAAAAwKdbupVOFzZw4EDp1KmTdO7c2cy7nZmZaaqZq9jYWKlbt64Zl20rvvbNN9/Yv//+++9l7969UrFiRWnSpIlH3wsAAAAAILB5Xeju3bu3nD9/XuLi4kzxtHbt2klSUpK9uFpKSoqpaG5z+vRpad++vX199uzZZvn5z38umzdv9sh7AAAAAADAK0O3GjFihFlcKRiko6KiJD8/301nBgAAAACAj47pBgAAAADAnxC6AQAAAACwCKEbAAAAAACLELoBAAAAALAIoRsAAAAAAIsQugEAAAAAsAihGwAAAAAAixC6AQAAAACwCKEbAAAAAACLELoBAAAAALAIoRsAAAAAAIsQugEAAAAAsAihGwAAAAAAixC6AQAAAACwCKEbAAAAAACLELoBAAAAALAIoRsAAAAAAIsQugEAAAAAsEgZqw4MAIDKzc2Vffv2yYULF6R69erSpk0bCQkJ8fRpAQAAuAWhGwBgmS1btsif/vQnOXv2rH1brVq1ZPjw4XLfffd59NwAAADcge7lAADLAveECROkUaNGkpCQIOvXrzdfdV236+MAAAD+jtANALCkS7m2cEdHR8vUqVOldevWUr58efNV13X7/PnzzX4AAAD+jO7lAIBSp2O4tUv5a6+9JsHBzvd3db1fv37y7LPPmv3at2/vsfME3CkrK0tSUlLc9nrJyclOX92lfv36Eh4e7tbXBABvRugGAJQ6LZqmGjZs6PJx23bbfkAg0MA9dOhQt79ufHy8W19v4cKF0qxZM7e+JgB4M0I3AKDUaZVydfz4cdOlvCDd7rgfEAi0BVgDaSC8TwDA/yF0AwBKnU4LplXKV65cacZwO3Yxz8vLM9tr165t9oP30YJ3s2bNMkME2rZtK/PmzZPOnTu73HfNmjUybdo0OXLkiFy9elWaNm0qL730kgwYMMC+T1BQkMvnvv766/Lyyy+b76Oiogp1g54+fbqMGTNG/IV2uaYFGAACD4XUAAClTufh1mnBtm/fLuPHj5cDBw7ITz/9ZL7qum4fNmwY83V7ocTERBk1apSpML9nzx4TumNiYuTcuXMu99feCuPGjTM/Ux2jP2jQILN8+umn9n3OnDnjtCxZssQE8V69ejkda/LkyU77Pffcc5a/XwAArBaUn5+fLwEuIyNDqlSpIunp6VK5cmVPnw4A+PU83drCrYE70Ofp9tZrT5cuXeRnP/uZvPPOO/aeCZGRkSYAF7fVuUOHDvLII4/IlClTXD7es2dPuXTpkmzcuNG+TVu6X3jhBbP4078nAMB/FffaQ/dyAIBlNFh37drVtIBq0TRtFdUu5bRwe6ecnBzZvXu3jB071r5NhwZ0797dtGTfiN7H37Rpkxw+fFhmzpzpcp/U1FRZt26dLFu2rNBjM2bMMEFdxwQ/+eST8uKLL0qZMq4/qmRnZ5vF8YMPAADeiNANALCUBmymBfMNaWlpZu70mjVrOm3X9UOHDhX5PL3DX7duXROC9eetvRsefPBBl/tq2K5UqZL85je/cdo+cuRI00KuN2a2bdtmgr92MZ8zZ47L4+h470mTJt3U+wQAwJ0I3QAA4JZoiN67d69cvnzZdBnXMeGNGjWS+++/v9C+Op5b52kvOI+zPsdGe0OEhobKM888Y8J1WFhYoeNoKHd8jrZ0azd4AAC8DaEbAAAYERERpqVau4A70nWtRl8U7YLepEkT8327du3k4MGDJiwXDN1ffPGF6XquxdqKM7b82rVrcuLECWnevHmhxzWIuwrjAAB4G6qXAwAAQ1uXO3bs6FTgTAup6Xp0dHSxj6PPcRxvbbN48WJzfK2IfiPacq5hvkaNGiV4BwAAeB9augEAgJ122R44cKB06tTJzM09d+5cyczMNNOAqdjYWDN+W1uylX7VfRs3bmyC9vr16+X999+X+fPnOx1Xu39/8MEH8sYbbxR6TS3StnPnTunWrZvpqq7rWkStf//+Uq1aNTe9cwAArEHoBgAAdr1795bz589LXFycmepNu4snJSXZi6ulpKSYFmgbDeQ6J/upU6ekXLly0qJFC1mxYoU5jqPVq1eb6uZ9+/Yt9JraTVwfnzhxognuDRs2NKHbccw2AAC+inm6mdsTAOABXHtKF/+eAABvvfYwphsAAAAAAIsQugEAAAAAsAihGwAAAAAAixC6AQAAAACwCKEbAAAAAACLELoBAAAAALAI83SLmHlDbSXfAQBwB9s1h5k7SwfXcgCAt17LCd0icunSJfM1MjLS06cCAAjAa5DO8Ylbw7UcAOCt1/KgfG6xS15enpw+fVoqVaokQUFBnj4dr7lrox9cTp48ed2J3hHY+D1BcfG7UphefvUiXadOHQkOZrTXreJaXhj/71Ac/J6guPhduflrOS3dOrA9OFjq1avn6dPwSvofiv9UuBF+T1Bc/K44o4W79HAtLxr/71Ac/J6guPhdKfm1nFvrAAAAAABYhNANAAAAAIBFCN1wKSwsTCZMmGC+AkXh9wTFxe8K4H78v0Nx8HuC4uJ35eZRSA0AAAAAAIvQ0g0AAAAAgEUI3QAAAAAAWITQDQAAAACARQjdfuapp56Snj17Fvn4119/Lb1795batWubIggNGjSQX/3qV/LJJ5+Yyd3ViRMnJCgoyL6EhoZKkyZNZOrUqfZ91MSJE83jPXr0KPQ6s2bNMo/df//9Fr1TuPv3yvb7ULZsWWnYsKG88sorkpWVZd/H8XfGttxzzz0ePW9Y9zdl8+bN5md88eJFs65/GxYtWiTR0dFm7s6KFStK69at5fnnn5cjR44U+rthW3Ruy3vvvVc+//xzp+NHRUXZ9ylfvrzceeed8t5777nhHQOex7UcVuBaHri4lnseoTuAfPTRR3LXXXfJ5cuXZdmyZXLw4EFJSkqSX//61zJ+/HhJT0932v9f//qXnDlzRr777juZNGmSxMfHy5IlS5z20Qv+Z599JqdOnXLarvvVr1/fLe8L7qEfyPT34dixY/Lmm2/Ku+++aypYOvrzn/9s9rEtH3/8scfOF+6jF+knn3xSRo4cKb/85S/ln//8p3zzzTeyePFiCQ8PNx/yHekF3PY7sn37dmnatKkJDAX/Bk2ePNns85///Ef69+8vQ4YMkX/84x9ufneAd+FajlvBtRxF4VpurTIWHx9eIjMzUwYPHiyPPPKIrFmzxumxli1bmscKFrK/7bbbpFatWuZ7vYuuf4T37Nlj9rWpUaOGdOzY0Vz4x40bZ7Zt27ZN0tLS5He/+535zwr/oK0ptt+HyMhI6d69u2zYsEFmzpxp36dq1ar2fRA4EhMTZfXq1SYMPPbYY/bt+mFdw0HBvy1lypSx/57oV70g69+Xb7/9Vn72s5/Z96tUqZJ9v1dffVVef/118zv38MMPu+29Ad6EazluFddyFIVrubVo6Q4Qerfqhx9+MN2IiqJdP4ry1Vdfye7du6VLly6FHvv9738vS5cudboz3q9fP9OVDf5J71bqBzJ+xlB/+ctfpHnz5k4X6eL+bcnOzjYXaf2Qp8dwJS8vT/7nf/5HfvzxR37nENC4lqM0cS2HI67l1iJ0Bwi966Qc/yN8+eWXZqyGbVm7dq3Tc+6++26zXf9j6B2rJ554QmJjYwsdW7uSZGRkyJYtW8xd+L/+9a/m4g3/or8f+vugXYx0TM65c+fk5Zdfdtqnb9++Tr9TH374ocfOF6X/s3dcHO9Q69+XghfZF154wb5vvXr1nB7bv3+//bFy5crJ7NmzzcVex4850jviuo+2zPz2t7+VatWqydNPP23xuwW8F9dy3Cqu5YGLa7ln0b08gLVp00b27t1rvtdxGNeuXSvUzUS7q129etXcDX3uuefMf5QZM2Y47afFOHSMht7h0jFCzZo1M8eGf+nWrZvMnz/ffBjTcWDarahXr15O++h27armOE4Q/vOzd7Rz507z/74o2kV1xIgRpgvstGnTnB7Ti7ptjOClS5fM3xrtwqpjSjt16mTfTz8IavEXHQum3w8fPtwUggLwf7iWoyS4lgcuruWeRegOEHohVocPHzbjMpTecbreL72O9bE9rhfso0ePymuvvWYqFuodUkd6N1y7q+kFnTvj/qlChQr23wftdti2bVtTXMNxXKCO2QnEP6SB9LO3cSy4pH9f9G+Lo9tvv90sOla0IFsVZZv27dublpS5c+fKihUr7NsjIiLMfrp88MEHplVGL+StWrUq5XcI+Aau5bhVXMsDF9dyz6J7eYB46KGHpHr16k6FMkoqJCTE3EHPyckp9JhWMNRFL9Ra+RD+LTg4WP74xz+aSrlXrlzx9OnAw7Qrol6otfjKrfx9ud7vkgYHnSJp7NixN/0agK/jWo7SxLUcjriWW4uWbj+kpfptXc0cq5fqvHj6i65VT3U6AL2jpVOO6FQjtv8ojrRYy9mzZ83FWcdtvPXWW6ZrSsGxGjabNm0y3de0iAL8n3Yh0m5CCQkJMnr0aE+fDjyoT58+puuZftULaUxMjNSsWVOSk5NNd7OCf1v0b4r+bXHskqbVkXXc1/XoPKF33HGHKQbl2HUN8Edcy+EOXMthw7XcWoRuP6QT3WsXD0fabUgv1FqlUu+QaxGVCxcumMns9RdepwjQIiqObON59D+ZjufROft0fs/rdVtB4NBxYDrOR6d+GDZsmKdPBx6kFU31Yrto0SIzHlR/J/RDuxZdeeCBB2TOnDlO+x84cMA+RrB8+fLSuHFjM87MVXEnR9oVTVv64uLiZP369Za+J8DTuJbDHbiWw4ZrubWC8gtOugYAAAAAAEoFY7oBAAAAALAIoRsAAAAAAIsQugEAAAAAsAihGwAAAAAAixC6AQAAAACwCKEbAAAAAACLELoBAAAAALAIoRsAAAAAAIsQugF4hRMnTkhQUJDs3bu32M+5//775YUXXrD0vAAAwI1xHQeKRugG/MBTTz1lLnR/+MMfCj327LPPmsd0H8f9e/bsWezjT5w40RyjR48ehR6bNWuWeUwvnAAAoOS4jgP+jdAN+InIyEhZvXq1XLlyxb4tKytLVq1aJfXr17/l49euXVs+++wzOXXqlNP2JUuWlMrxAQAIZFzHAf9F6Ab8RIcOHcwFe82aNfZt+r1eSNu3b3/Lx69Ro4Y89NBDsmzZMvu2bdu2SVpamjzyyCNO++bl5cnkyZOlXr16EhYWJu3atZOkpCSnfXbt2mXOKzw8XDp16iRff/11odf8z3/+Iw8//LBUrFhRatasKQMGDDCvBwCAv+E6DvgvQjfgR37/+9/Ln//8Z6e714MGDSrV4y9dutTp+P369ZPQ0FCn/d566y154403ZPbs2bJv3z6JiYmRxx57TL777jvz+OXLl+VXv/qVtGrVSnbv3m26vY0ePdrpGBcvXpRf/OIX5oL+1VdfmYt9amqqPPHEE6X2fgAA8CZcxwH/ROgG/Ej//v1l69atkpycbJZ///vfZltp0QtsRkaGbNmyRTIzM+Wvf/2ruYAXpBfpV199Vfr06SPNmzeXmTNnmrvkc+fONY9rVzm9i7548WJp3bq1Oe7LL7/sdIx33nnHXKinTZsmLVq0MN/rhwPtGvftt9+W2nsCAMBbcB0H/FMZT58AgNJz++23my5iehc7Pz/ffB8REVFqxy9btqy5+Otd+GPHjkmzZs2kTZs2Tvvoxfz06dPStWtXp+26/r//+7/m+4MHD5rnaZc0m+joaKf9dV+9MGuXtIKOHj1qXhsAAH/CdRzwT4RuwM/oHesRI0aY7xMSEiw5fpcuXcw4LVd3x0uLdl179NFHzd11V8VgAADwR1zHAf9D93LAz+h0IDk5OXL16lUzBqu0aTcyXfRi/eSTTxZ6vHLlylKnTh3TJc6RruvYL9WyZUszRkyrstrs2LGjUEGZAwcOSFRUlDRp0sRpqVChQqm/LwAAvAHXccD/ELoBPxMSEmK6fX3zzTfm+6Kkp6fL3r17nZaTJ08W6zU2bdokZ86ckapVq7p8XMd16Z3txMREOXz4sIwZM8Yc//nnnzeP60Ve5wQdMmSIOc/169eb8WMF5yW9cOGC9O3bV7788kvTFe3TTz81BWVyc3NL9G8CAICv4DoO+B+6lwN+SO9S38jmzZsLTUEyePBgee+992743BvdoR45cqT5MPDSSy/JuXPnzJ3xjz/+WJo2bWoe1/Fdn3zyifzhD38w56CP68W9V69e9mPY7rJrIRed4iQ7O1saNGhgWgCCg7lfCADwX1zHAf8SlK9VGgAAAAAAQKnjNhMAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAAGARQjcAAAAAABYhdAMAAAAAYBFCNwAAAAAAFiF0AwAAAABgEUI3AAAAAAAWIXQDAAAAACDW+H8RQvlBTbsN+gAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 1000x500 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure, axes = plt.subplots(1,2,figsize=(10,5))\n",
    "ax_0 = sns.boxplot(x=\"model\",y=\"R2\",data=model_combo_df,ax=axes[0])\n",
    "ax_0.set_xlabel(\"ML Model\")\n",
    "ax_0.set_ylabel(\"$R^2$\")\n",
    "ax_1 = sns.boxplot(x=\"model\",y=\"MAE\",data=model_combo_df,ax=axes[1])\n",
    "ax_1.set_xlabel(\"ML Model\")\n",
    "plt.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfd76d5f",
   "metadata": {},
   "source": [
    "Tukey HSD test was performed on the results from the three regressors. There is no significant difference for $R^2$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "6d81ca33-91c1-42ec-92ad-fd59739724e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multiple Comparison of Means - Tukey HSD, FWER=0.05\n",
      "===================================================\n",
      "group1 group2 meandiff p-adj   lower  upper  reject\n",
      "---------------------------------------------------\n",
      "  HGBR   LGBM  -0.0505 0.3259 -0.1345 0.0334  False\n",
      "  HGBR     RF   0.0331 0.6154 -0.0509 0.1171  False\n",
      "  LGBM     RF   0.0836 0.0513 -0.0004 0.1676  False\n",
      "---------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "model_tukey_r2 = pairwise_tukeyhsd(endog=model_combo_df[\"R2\"], groups=model_combo_df[\"model\"], alpha=0.05)\n",
    "print(model_tukey_r2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "41ec5304",
   "metadata": {},
   "source": [
    "Tukey HSD test was performed on the results from the three regressors. There is no significant difference for MAE. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "19e93140-fa5c-43f8-9693-df12c6e2f849",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Multiple Comparison of Means - Tukey HSD, FWER=0.05\n",
      "===================================================\n",
      "group1 group2 meandiff p-adj   lower  upper  reject\n",
      "---------------------------------------------------\n",
      "  HGBR   LGBM   0.0065 0.8527 -0.0225 0.0356  False\n",
      "  HGBR     RF  -0.0112 0.6286 -0.0402 0.0179  False\n",
      "  LGBM     RF  -0.0177 0.3162 -0.0468 0.0113  False\n",
      "---------------------------------------------------\n"
     ]
    }
   ],
   "source": [
    "model_tukey_mae = pairwise_tukeyhsd(endog=model_combo_df[\"MAE\"], groups=model_combo_df[\"model\"], alpha=0.05)\n",
    "print(model_tukey_mae)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8765f74b",
   "metadata": {},
   "source": [
    "Determine the model with the highest $R^2$ and lowest MAE. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "1a16cbf4-0d05-4b76-a309-deb6b92f0c94",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('RF', 'RF')"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_best_r2 = model_combo_df.groupby(\"model\").mean().reset_index().sort_values(\"R2\",ascending=False).model.values[0]\n",
    "model_best_mae = model_combo_df.groupby(\"model\").mean().reset_index().sort_values(\"MAE\",ascending=True).model.values[0]\n",
    "model_best_r2,model_best_mae"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2a80860",
   "metadata": {},
   "source": [
    "Evaluate the statistical significance of the $R^2$ and MAE comparisons across the three models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "bbe58023-9e11-4976-a3bb-75cc3ab6573f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA00AAAIoCAYAAAC1TQBxAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAQ8NJREFUeJzt3QmYHFW5OO4zTEgIhJkYtiQQQghbQNawiSI7xg1wQQRlUcQrV8V7uYqCIuByA6I/8SqisopcNpXlihqFyyKIIIvsXDbZJYAEZhKWhEz6/3zlv9uZyUwlaVLp6er3fZ5Karpruk9Vf3NOf1WnzmmrVCqVBAAAwICWG/hhAAAAgqQJAAAgh6QJAAAgh6QJAAAgh6QJAAAgh6QJAAAgh6QJAAAgh6QJAAAgh6QJAAAgh6QJAAAgh6SJhvnhD3+Y2tra0nbbbTfoNvF8LN/5zncWeu6cc87Jnrv11lsXeu6Pf/xjet/73pfWWGONNGLEiLTOOuukf/mXf0lPPPFEbZvnnnsuDRs2LH30ox8d9P1nz56dRo4cmd7//vcvcdkb6fbbb0977bVXGjNmTFpxxRXTm9/85vRf//Vfi/y9Sy65JO23335p3XXXzX5vww03TP/xH/+RXnrppUGPz1FHHZUmTZqUHec111wzffCDH0yvvPLKUt2fuXPnpi9+8Ytp/Pjx2ecRx/3KK6/ss02856mnnpr23HPPNG7cuLTyyiunLbfcMp122mmpp6dnqZYHKC9t09J37733pn333bfWtqy66qrp7W9/e/rVr35V1+sddthh2X6+5z3v6fP4Cy+8kE4++eTstVdbbbU0evTotP3226eLLrooFWHBggXpW9/6VtYGrrDCCmmzzTZLF1xwwULb/fnPf07/+q//mqZOnZqWX375rOw0oQo0yA477FBZZ511KhGGDz300IDbxHOxrLHGGpWXX365z3Nnn3129twtt9zS5/H/+q//qrS1tVUmT55c+frXv14544wzKv/xH/9R6ezszJY//vGPtW2nTZtWGTVq1EKvXXXOOedk7/HLX/5yicveKL/73e8qw4cPr2y33XaV//f//l/lJz/5SeWLX/xi5Qtf+MIif3eVVVapbLrpppVjjz22cvrpp1eOOOKI7LU22mijyiuvvNJn25deeqmy+eabZ79z9NFHV84888zKiSeeWHn3u99dmTVr1lLdpw9/+MOVYcOGVT7/+c9XfvzjH1fe8pa3ZD9ff/31tW3uvvvu7HPffffdK9/61rcqP/rRjyrve9/7ss/ooIMOWqrlAcpL27T0/frXv6684x3vqBx//PFZm3TKKadUdtxxx6ycUacviTiuUf+vsMIKWXvT269+9avK8ssvX9l7772z9/jBD35Q2WWXXbL3+epXv7qU96pS+dKXvpS99mGHHZbtV5Qnfr7gggv6bHfcccdl5Zo6dWplgw02yLah+fjUaIi//vWvWaVxySWXVFZbbbWsIh1IbLPFFltk/3/nO99ZZMN0ww03VJZbbrmsMu7f2Dz88MNZAzdu3Ljal/qf/exnA1ZwVXvuuWfWmL322mtLXPZG6OrqyvYxkoWenp4l/v1rrrlmocd++tOfZvsbSVRvhx9+eGX06NHZ8SjSzTffnL3/ySefXHvs1Vdfzb54RPJU9fzzz1fuueeehX7/Yx/72JD7AgEMTdqmZWf+/PnZibcNN9xwsX9nwYIFWb3/8Y9/vDJx4sSFkqY4Bo899thCv7PrrrtWRowYUZkzZ85SK/9TTz2VJUKf/vSn+7xXfMZrrbVWtn9VM2fOrJ14jO0lTc3Jp0ZDxFm2N73pTZW5c+dmX77XX3/9AbeLiiUqmKjwolHpfbVjoIYpzmS1t7cP+kW+mgBMnz49+zkq0JVWWqny3ve+d6Ftn3322ey1Dj300LrKPpBqJR9Xg6KxiEp8ypQpC50trNdpp52W7d99991X2796kqfeuru7s9c88sgja4+9+OKL2Vm+o446Kvs5jkXvxru/+++/v/KBD3wgO26xz3G27fLLL1+s948rZPE5RELY23/+539m5XriiSdyf/9//ud/su3if4A82qZi2qbBvOc978mO3+KK47TyyitXnnnmmQGTpsHEVb44vnfddddCiU+cWFt99dWzXhUbb7xx1mticZx66qnZa9577719Hj///POzx3v3hOhN0tS83NNEQ/z3f/931hd7+PDhaf/9908PPfRQuuWWWwbd/vjjj0/PPvtsdn/KYOKelv/93/9NO+64Y9a/eCBxv070I7/iiiuyn1daaaW09957p9/97ndp1qxZfbaNPtBxL8xHPvKRN1T2/mL7KMc73/nONH369KzvevT17n2PTvST/vvf/75Yy+uvv177vauuuip1dHSkp59+OrsfadSoUdnPhx9+eHrttddSPWbOnJn9H33Qq2644Ybs9dZbb73sHqboox7969/61remO+64Y6G+7NGn/P77709f+tKXsnsA4rjvs88+6dJLL13k+//lL39JG2ywQbYfvW277bbZ//3fb3HKDzAQbVMxbVPVyy+/nD33yCOPpO9+97vpt7/9bdptt90Wq3xxH1fc23rMMceksWPHpiUxUDsQn1u0TdFufuYzn0nf+973sjbt0EMPTaeccspitU3xOU2ZMmXAtimep2QanbXRem699dbsLMuVV15Zu5wdl7I/97nPDXo2L0S/5LFjx9bO6PU/m3fHHXdkPw/0Or1tttlmlTFjxvTpaz1Qv+rtt9++suaaa/a5UrMkZR9InBnr3w89rqBEt4wtt9yy9tijjz5a6zO/qKV3l7rYtxVXXDFbPvvZz2bvE//HdnFfUD3ibGac1XzwwQdrj8W9UvGacT/TtttuW/nv//7vyg9/+MPsjGGc6fzb3/5W23a33XbL7pPqfSUqjlv0vV+cM6GbbLJJdja3vzi7F2WIe5cGE2dc48zhpEmTKq+//voS7jnQSrRNxbVNVf/yL/9Sez66K37wgx9c7Htg457WqMurbcniXml64YUXsitJ0W2uf9sW+/f3v/+9z+PRVkbXx/738fYX773uuusu9Hh0v4z9i/udBuJKU/Ma1uikjdYTZ8Ni5KBddtkl+zlGkYmzW+edd152FaK9vX3QM3o77bRT+tGPfpT+/d//fcCzUCFGTcsTz3d3d9d+jtHWYpSd888/P33yk5/MHnv00UfTTTfdlD7/+c+n5ZZb7g2XvbcYAS5GT6qKKygHHXRQOumkk7KzYXEGLZb+o8MNZvPNN6+tz5kzJzur+alPfao2Wl6ceZw3b1768Y9/nL72ta+l9ddfPy2uOCZnnnlmNkJe79+L96nuf5xBjStaIUare8tb3pKNYveNb3wjO0N69dVXZ+8bn0/1MwrveMc70nHHHZddFYtR9wbz6quvZmdg+4uRiqrPDybOHt53333p17/+dXbWFGAw2qbi2qaqf/u3f8t6J/ztb39LF198cXbFLNqnRXnwwQezK0ExMt1A7cFg4spYXJGLEWC///3v1x6PvPeXv/xl+tCHPpStx9Wv3m3ThRdemI1CG70nimibaE6+RbBMRQUZlVFU7FH5V8XwqFGxxxfwaCgGEkOIxu/F8J6RFPRXbZB6fzEfSDzfu/GKL9PRuMRQrdUv8NFIhd7dH95I2XuLy//9hxuN7mfhscceyxqlqHR33333tKSii1yIrhm9HXDAAVnS9Kc//Wmxk6brr78+66YQDcg3v/nNAd/nve99by1hCtHVIbqf3HjjjdnPDz/8cNYgHXvssdkykBheN/b5+eef7/N4DJce3UzivWLI8f6q3Q2rZekvhp09/fTT09e//vX0rne9a7H2GWhN2qZi26aqjTbaKFtCJGRRrmhHbr755txhuD/3uc+lHXbYIX3gAx9Yovf77Gc/m2bMmJHOPffcPklctDeRSP3kJz/JlsHapt5d+6o6OzuzdqfetonmJWlimYqrDs8880xWwcfSX5wty6vc48rEzjvvnCUAMf9C/wo/Gpm77rpr0N+PCu6BBx5IW2+9dZ/HYz6MH/zgB9lZrDiDF/9vvPHGaYsttlhqZV8S0Qj2TyIGU00uqmcK4x6iOOPY2+qrr579/+KLLy7Wa955553ZPE8xv9MvfvGLha7SxPuE/u9Tfa/q+8RZvhDHNJKvgcTn9uSTTy7U1/+aa67JPuuYcym+MPQXn0XvsvSfJyX6vscXmK985SuLtc9A69I2Fds2DSauOsU8VXElKe7DHUjsXyQ+MY9gJG9V8+fPz67mxGPxXv3vez3hhBOyhPPEE09MBx54YJ/nqm1THN+DDz54wPeNOZdCtEG9nX322emQQw7JHo92Kk4M9k748tommpukiWUqKu/4Uh3dt/qLCjEGBoguDoOdoYkuENEwRXeBr371q32eixsy40xbVLCPP/54mjhx4kK/H90BonHqPyFenJGbPHlydhZvjz32yBKP/ldX3mjZq6pXX3pXstFghJjoMAyURAymmlyEmDgvuk5UB4Koiq4QIbp6LErcoDtt2rRsX3/zm9/0uZJUFe8TBkpm4r2qZxJjIsMQk/nlnZ2M5/t3+aieFYwvB7GP0W2ld6MYZyarz/d2+eWXp0984hNZt8SBPiuA/rRNxbZNg6l2Yevq6hp0m+rEv/0n8q22QVGeGFQiuv5VxbGIbpPxWJxA6y/awriqF0ngoq6c9W+bNtlkk1rbc8YZZ2SDHEUiu6i2iRJo9E1VtI64qTKGCo35FQYSE/tFSF544YUD3mxbde211/aZI6P3sK7XXXdddnPpzjvvvNBNnDHUa9ys23sujN5i4rt4vZgULyYgjBte30jZl/Rm29if3vMQxQ29i7P03pfbb789e/0DDjigz/vuv//+2WSATz/9dO2xxx9/PBsKvLcYxjVubB0/fnyf/R9IDEvb0dGRzY9UFcPVxvvH5LJV8VnEzc29B4eoeu655yqLctNNNy00T1PcCLzeeutlE/j2Fp9/DIUeN2bnDYEOUKVtKr5timHS+5s3b15lq622qowcObIye/bs2uPRVkTbFM9X26pLL710oSXmotp6662z9Zjrqir2NY71Rz7ykWxAjMEccsgh2TDjMTF6PW3Tk08+Oeg8TTFQR+95mnozEETz8qmxzERFFhXFZZddNuDzMRJQVIK956UYqGEKO+20U20Env6zrldHdosv1d/4xjeyORdirp+YiDW+5MckgwOJ0eGqr/nWt771DZd9sIYpZgOPssTIOt/97nezkeWigp8xY0ZlaYjGM8r6oQ99KJtHYt99981+Pvroowc8hv0ToXgs5l+KyRV7L7///e/7bHv11Vdno+rFxIRxzGPG82i8Y/96N4Axyl2MqBcj7cU+x6zpMZ/Iu971rmy0qMUR+xBJX3yOMZJUjLwXP8cXkaqY0DBGPIoGOPa7f/nvvPPOOo8oUGbapuLbpn322ScbBTUm3I2J0qMN2GijjQacHPjggw/OHl/UibuBRs+LydAjEYp9PuussxZqBx555JE+E87Ga8RoszHKYLQtMU9WtDfRZi2O+PyirJ/85Cez/YryxM8xomxv0T7FPscSJ/tim+rP55577mK9F40naWKZiUo7rgL0nw29/5mfOHNTHQJ0sIYphjIdrGEKf/jDH7Kzcquuumr2emuvvXblsMMOW2im8P622Wab7DVj+Ow3WvZFTSAYCUNMIBgNx89//vPK0hJn56JhiveK8kQDHQ1gfwMlTXnDx8b2/cXZxBj+No5NXE068MADs6tV/UVDddBBB2VnU6NMcRYuJjX8xS9+sVj7FGc3Y7jZ+P04ZvE59W/Ie8fEQEskdQD9aZuKb5suuOCCyu67755NSxEnvCIpiZ8HmuT8jSRN1eHeB1vi+d7iClh8jhMmTMiOUbQxMU1GnNxbHJGUxkTrUZZI1mKKjPPOO2+h7fLap4HaVoamtvin0V0EoVVEv/AYXKE6gSEANJq2CRbtn4P8AwAAsBBJEwAAQA5JEwAAQA73NAEAAORwpQkAACCHpAkAACDHsFRic+fOzZaqBQsWpFmzZqVVVlkltbW1NbRsAK0keoLPnj07jR8/Pi23XGufr9M2ATRf21TqpGn69OnphBNOaHQxAPj/Pfnkk2mttdZKrUzbBNB8bVOpB4Lofzavq6srrb322tmB6ejoaGjZYEnEmeiI2zBhwoSWP1NP88VSd3d39n4vvfRS6uzsTK2s0W2T+oSyEtsU2TaV+krTiBEjsqW/aJQkTTST+fPnp8cffzxbnzJlSho2rNR/upQ4lnQ/a3zb1OgYgKKIbYpsm6TgAAAAOSRNAAAAOSRNAAAAOSRNAAAAOSRNAAAAOSRNAAAAOYzFCE0g5prYZpttautQL7GEGKCsxDZFkjRBE4jKf4011mh0MSgBsYQYoKzENkWShgMAAORwpQmawIIFC9LTTz+dra+55pq6HVA3sYQYoKzENkWSNEGTNAR33nlntj5u3DgNAXUTS4gBykpsUyTRBAAAkEPSBAAAkEPSBAAAkEPSBAAAkEPSBAAAkEPSBAAAkMOQ49AEYtjUrbbaqrYO9RJLiAHKSmxTJEkTNIGo/MePH9/oYlACYgkxQFmJbYokDQcAAMjhShM0ySznM2fOzNbHjh2r2wF1E0uIAcpKbFMk0QRN0hDcfvvt2RLrUC+xhBigrMQ2RZI0AQAA5JA0AQAA5JA0AQAA5JA0AQAA5JA0AQAA5JA0AQAA5DBPEzSBmGti8803r61DvcQSYoCyEtsUSdIETSAq/wkTJjS6GJSAWEIMUFZimyJJwwEAAHK40gRNIGY2f/7557P11VZbTbcD6iaWEAOUldimSKIJmqQhuOWWW7Il1qFeYgkxQFmJbYokaQIAAMghaQIAAMghaQIAAMghaQIAAMghaQIAAMghaQIAAMhhniZoAjHXxJvf/ObaOtRLLCEGKCuxTZEkTdAEovJfZ511Gl0MSkAsIQYoK7FNkaThAAAAOVxpgiZQqVTSCy+8kK2vssoqqa2trdFFokmJJcQAZSW2KZIrTdAEenp60k033ZQtsQ71EkuIAcpKbFMkSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAO8zRBk8xyPmXKlNo61EssIQYoK7FNkSRN0ASi8p88eXKji0EJiCXEAGUltimSNBwAACCHK03QBCqVSurq6srWOzs7U1tbW6OLRJMSS4gBykpsUyRXmqAJ9PT0pBtuuCFbYh3qJZYQA5SV2KZIkiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAc5mmCJpnlfP3116+tQ73EEmKAshLbFEnSBE0gKv8NN9yw0cWgBMQSYoCyEtsUSRoOAACQw5UmaAKVSiXNmTMnWx81alRqa2trdJFoUmIJMUBZiW2K5EoTNIGenp503XXXZUusQ73EEmKAshLbFEnSBAAAkEPSBAAAkEPSBAAAkEPSBAAAkEPSBAAAkEPSBAAAkMM8TdAks5yvu+66tXWol1hCDFBWYpsiSZqgCUTlv/HGGze6GJSAWEIMUFZimyJJwwEAAHK40gRNoFKppFdffTVbHzlyZGpra2t0kWhSYgkxQFmJbYrkShM0gZ6ennT11VdnS6xDvcQSYoCyEtsUSdIEAACQQ9IEAACQQ9IEAACQQ9IEAACQQ9IEAACQQ9IEAACQwzxN0ARiromJEyfW1qFeYgkxQFmJbYokaYIm0N7enjbddNNGF4MSEEuIAcpKbFMk3fMAAAByuNIETaBSqaR58+Zl68OHD9ftgLqJJcQAZSW2KZIrTdAEenp60pVXXpktsQ71EkuIAcpKbFMkSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAO8zRBE4i5JtZaa63aOtRLLCEGKCuxTZEkTdAE2tvb0xZbbNHoYlACYgkxQFmJbYqkex4AAEAOV5qgCVQqldrs5nEmTbcD6iWWEAOUldimSK40QROIRmDGjBnZUm0QoB5iCTFAWYltiiRpAgAAyCFpAgAAyCFpAgAAyCFpAgAAyCFpAgAAyCFpAgAAyGGeJmgCMdfEuHHjautQL7GEGKCsxDZFkjRBE4hJ+qZOndroYlACYgkxQFmJbYqkex4AAEAzJE2HHHJIdik1luWXXz5NmjQpHXXUUem1116rbVN9vvfytre9raHlBoaGmPz92mtTuuCCf/xvMnig0dRLUB5DqnvetGnT0tlnn51ef/31dNttt6WDDz44S4xOOumk2jbxfGxXNXz48AaVFpad+fPnpxkzZmTrEf/Dhg2pP92Gu+SSlD73uZSeeuqfj621Vkrf+15K739/I0s29IglxMCyoV5a9sQ2LXGlKYwYMSKNHTs2TZgwIe2zzz5p9913T1deeWWfbUaPHp1tU13GjBnTsPICQ+OLyQc/2PeLSXj66X88Hs8DLEvqJSifIZuC33PPPenGG29MEydObHRRYEh55ZW42bXRpRgaoqvLEUekVKks/Fw8FoMnxZne3Xd3zKp0D6I39cnSp15qHPUbLZM0XXHFFWnUqFHZ5dW5c+em5ZZbLv3gBz/os83++++fjY5Sdd5552VXpQYSrxFLVXd3d4Glh2Vj9dUjthtdiuYQX1DiTG9nZ6NLMnSMGJHSL37R6FK0tqHUNqlPlj31UnHUb7RM0rTLLruk0047Lb388svpu9/9btYX9QMf+ECfbeLx6LZXVR2PfyDTp09PJ5xwQqFlBoAloW0CaD5DKmlaaaWV0nrrrZetn3XWWWnzzTdPZ555Zjr00ENr28R9TNVtFuXoo49ORx55ZJ+zeXG/FDSz557TpaPqD39I6V3vWvR2v/lNSm9/+7IoUXN0X4njRuMMpbZJfbL0qZcaR/1GyyRNvUXXvGOOOSZrWA444IA0cuTIugaWiAXKZMUVUzIg0D/suec/RqOKm6sHun8g7h2I52M7Xwz/Yf78RpeAodQ2qU+WPvVS46jfaJnR8/rbd999s/uXTj311EYXBRoqht5fffXVsyXW+Yf4whHD94b+h6X68ymn+GLSm1hCDBRLvdQ4YpuWTZrinqbPfOYz6Vvf+lZ2nxO0qjh5sO2222ZL74FQ+Md8J3Hj75pr9n08zuTG4+ZD6UssIQaKp15qDLFNkdoqlYEuHpdT9Bvv7OxMXV1dqaOjo9HFAZZyX/brr0/pmWdigJiUdtzRmdyhRP07OMemvNRLUJ76V09moBTii8jOOze6FAD/pF6C8pA0QROIucuuvPLKbH2PPfbIuq5CPcQSYoCyEtsUSTRBk+gx1TlLiVhCDFBWYpuWHAgCAACg0SRNAAAAOSRNAAAAOSRNAAAAOSRNAAAAOYyeB02gra0tjRkzprYO9RJLiAHKSmxTpLZKpVJJLcKs6wCNof4dnGMDMPTrX93zAAAAckiaAAAAcrinCZrA/Pnz09VXX52t77rrrmnYMH+61EcsIQYoK7FNkUQTNIl58+Y1ugiUhFhCDFBWYpui6J4HAACQQ9IEAACQQ9IEAACQQ9IEAACQQ9IEAACQw+h50ATa2tqyGaur61AvsYQYoKzENkVqq1QqldQiuru7sz+mrq6u1NHR0ejiALQM9e/gHBuAoV//6p4HAACQQ9IEAACQwz1N0AR6enrStddem63vvPPOqb29vdFFokmJJcQAZSW2KZKkCZpA3Hr46quv1tahXmIJMUBZiW2KpHseAABADkkTAABADkkTAABADkkTAABADkkTAABADqPnQRNoa2tLo0aNqq1DvcQSYoCyEtsUqa3SQmMydnd3p87OztTV1ZU6OjoaXRyAlqH+HZxjAzD061/d8wAAAHJImgAAAHK4pwmaQE9PT7r++uuz9R133DG1t7c3ukg0KbGEGKCsxDZFkjRBE4hbD+fMmVNbh3qJJcQAZSW2KZLueQAAADkkTQAAADkkTQAAADkkTQAAADkkTQAAADmMngdNoK2tLY0cObK2DvUSS4gBykpsU6S2SguNydjd3Z06OztTV1dX6ujoaHRxAFqG+ndwjg3A0K9/dc8DAADIIWkCAADI4Z4maAI9PT3pxhtvzNZ32GGH1N7e3ugi0aTEEmKAshLbFEnSBE0gbj2M/rbVdaiXWEIMUFZimyLpngcAAJBD0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJDD6HnQJIYPH97oIlASYgkxQFmJbYrSVmmhMRm7u7tTZ2dnNhxlR0dHo4sD0DLUv4NzbACGfv2rex4AAEAOSRMAAEAO9zRBE+jp6Uk333xztr7ddtul9vb2RheJJiWWEAOUldimSJImaAJx6+GsWbNq61AvsYQYoKzENkXSPQ8AACCHpAkAACCHpAkAACCHpAkAACCHpAkAACCH0fOgSRg6laVFLCEGKCuxTVHaKi00JmN3d3fq7OxMXV1dqaOjo9HFAWgZ6t/BOTYAQ7/+1T0PAAAgh6QJAAAgh3uaoAn09PSk2267LVufOnWqPtvUTSwhBigrsU2RJE3QBOLWw+eee662DvUSS4gBykpsUyTd8wAAAHJImgAAAHJImgAAAHJImgAAAHJImgAAAHJImgAAAHK0VVpoTMbu7u7U2dmZurq6UkdHR6OLA9Ay1L+Dc2wAhn7960oTAABADkkTAABAjmF5TwJDQ09PT7rjjjuy9S222CK1t7c3ukg0KbGEGKCsxDZFcqUJmkDcevjMM89kSwvdhkgBxBJigLIS2xRJ0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJCjrdJCYzJ2d3enzs7O1NXVlTo6OhpdHFhs8Wca80+EmHeira2t0UWiSTUqltS/Q+fYqE8oK7FNkfWvyW2hCUTFP2yYP1feOLGEGKCsxDZF0j0PAAAgh3QcmkB0N7j77ruz9U033TTrdgD1EEuIAcpKbFMkV5qgSfppP/XUU9nSQrchUgCxhBigrMQ2RZI0AQAA5JA0AQAA5JA0AQAA5JA0AQAA5JA0AQAA5JA0AQAA5GirtNCYjN3d3amzszN1dXWljo6ORhcHFlv8mc6bNy9bHz58eDbrOTRTLKl/h86xUZ9QVmKbIutfk9tCE4iKf8SIEY0uBiUglhADlJXYpki65wEAAORwpQmaQE9PT7rvvvuy9Y033ji1t7c3ukg0KbGEGKCsxDZFcqUJmqSf9uOPP54tLXQbIgUQS4gBykpsUyRJEwAAQA5JEwAAQA73NNFwCxYsSE888USaPXt2WnnlldPaa6+dlltOPg8ANI7vJ/S2xJ/8IYcckvbZZ59Bn//LX/6S9ttvvzRu3Lhs2MeJEyem97znPelXv/pVrX/pY489lg0LWV1iLP311lsvfeMb3+jTB/X444/Pnp82bdpC73PyySdnz+28885LugsMIffff3/63ve+l37605+mSy65JPs/fo7HAQAawfcT+luq6fLll1+ett9++zRnzpwsuCKwZsyYkd73vvelr3zlK9nEUb1dddVV6ZlnnkkPPfRQOuGEE9I3v/nNdNZZZ/XZJpKva665Jj311FN9Ho/tIuOneUV8XHzxxdnEYr3Fz/G4igkAWNZ8P6HQ7nkvv/xyOvTQQ9O73/3uLCPvbcqUKdlz/UcyWWWVVdLYsWOz9bgidfbZZ6fbb78927Zq9dVXT1OnTs2SsC9/+cvZYzfeeGP6+9//nvbdd9/a0JJDWXV2avpe8v7tb3+bu008P2nSJJfCU0rz58/PhlIFAN8rFt1WxjGK7xpLyveTxho+fHgqfdL0+9//Pr3wwgvpqKOOGnSb6E43mFtvvTXddttt6aCDDlrouY9//OPZ61aTprjK9JGPfGSRZZo7d262VPU/Y7CsTJ8+vSHv2+yiD/FJJ53U6GIMOXvttVeji0ATi3lLdt1119o6y16j2yYx0Px8r1i0O+64o7DX9v2kOMcdd1waqpZaivzggw9m/2+44Ya1x2655ZY0atSo2nLFFVf0+Z0ddtghezyyym222SZ96EMfGjBpinuiolH5wx/+kF3RikujkUgtTqXS2dlZWyZMmLBU9hUaKe/kAyxO/Ky44orZIpYao9FtkxgAGGKj52222Wa1TH/99dfPLpv2dtFFF2Vd915//fV0zz33pM9+9rPpTW96UzrxxBP7bLf88sunj370o1n3vb/+9a9pgw02yF57UY4++uh05JFH1n6OxKsRiVOUg75i4rnzzz9/kdsdcMABWddNgLIYKm0Tzcv3iuL4fkLhSVMkReGBBx7IBoMIMXpejIo3mGgkqs9H8vTII4+kY489Nhs1b4UVVuizbVxZ2m677bLkanGuMlXfP5ZGG8r9Mxtl8uTJqaOjI7dbSjwf2+kz/I8+1v/3f/+XrW+00UaOCXUTS43X6LZJDDQ/3yuKi23fTxjMUvu099xzzzRmzJg31Mcz+lbH1aiBbnDcZJNNsiWSpsjuaW5R0Qw0lHxv8bwK6Z8NQVxljaWeG1uhSiwhBiirpRHbvp+wVK80xdDh/W+wi5HwzjjjjGyOphhB74gjjsiuPsXw4zHs+EA3nMbAETNnzswSpbvvvjsb/36XXXbJMviBXH311VlXvtGjR9dTbIaYuLoY97FFfPQ+oxOff1RI8TwAwLLk+wlLLWm69tpr05ZbbtnnsRgmPJKmGA48rjbFgA6zZs3KbnLdeuut04UXXpgN6NDb7rvvXkumYj6md73rXdlcTYNZaaWV6ikuQ1hUPDF4iBm3AYChwvcT3nDSdM4552TLYCJB+vnPf577Guuss85CczYNJO5timUwp5xyyiJfg6EvKqCICQCAocL3E3qTLgMAAOSQNAEAAOSQNAEAADRqcltg6YjBUnbaaafaOtRLLCEGKCuxTZEkTdAE2traspF74I0SS4gBykpsUyTd8wAAAHK40gRNIGY2f+ihh7L1mDTaPBHUSywhBigrsU2RJE3QZA3B5MmTNQTUTSwhBigrsU2RRBMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOQ45DE2hvb09ve9vbautQL7GEGKCsxDZFkjRBE2hra0ujR49udDEoAbGEGKCsxDZF0j0PAAAghytN0CSznD/66KPZ+qRJk8xyTt3EEmKAshLbFEnSBE3SENx///3Z+sSJEzUE1E0sIQYoK7FNkUQTAABADkkTAABADkkTAABADkkTAABADkkTAABADkkTAABADkOOQxNob29P22+/fW0d6iWWEAOUldimSJImaAJtbW1p1VVXbXQxKAGxhBigrMQ2RdI9DwAAIIcrTdAks5w/8cQT2fraa69tlnPqJpYQA5SV2KZIkiZokobgnnvuydbXWmstDQF1E0uIAcpKbFMk0QQAAJBD0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJDDkOPQBGLY1G222aa2DvUSS4gBykpsUyRJEzSBqPzXWGONRheDEhBLiAHKSmxTJGk4AABADleaoElmOX/66aez9TXXXFO3A+omlhADlJXYpkiSJmiShuDOO+/M1seNG6choG5iCTFAWYltiiSaAAAAckiaAAAAckiaAAAAckiaAAAAckiaAAAAckiaAAAAchhyHJpADJu61VZb1dahXmIJMUBZiW2KJGmCJhCV//jx4xtdDEpALCEGKCuxTZGk4QAAADlcaYImmeV85syZ2frYsWN1O6BuYgkxQFmJbYokmqBJGoLbb789W2Id6iWWEAOUldimSJImAACAHJImAACAHJImAACAHJImAACAHJImAACAHJImAACAHOZpgiYQc01svvnmtXWol1hCDFBWYpsiSZqgCUTlP2HChEYXgxIQS4gBykpsUyRpOAAAQA5XmqAJxMzmzz//fLa+2mqr6XZA3cQSYoCyEtsUSTRBkzQEt9xyS7bEOtRLLCEGKCuxTZEkTQAAADkkTQAAADkkTQAAADkkTQAAADkkTQAAADkkTQAAADnM0wRNIOaaePOb31xbh3qJJcQAZSW2KZKkCZpAVP7rrLNOo4tBCYglxABlJbYpkjQcAAAghytN0AQqlUp64YUXsvVVVlkltbW1NbpINCmxhBigrMQ2RXKlCZpAT09Puummm7Il1qFeYgkxQFmJbYokaQIAAMghaQIAAMghaQIAAMghaQIAAMghaQIAAMghaQIAAMhhniZoklnOp0yZUluHeoklxABlJbYpkqQJmkBU/pMnT250MSgBsYQYoKzENkWShgMAAORwpQmaQKVSSV1dXdl6Z2dnamtra3SRaFJiCTFAWYltiuRKEzSBnp6edMMNN2RLrEO9xBJigLIS2xRJ0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJBD0gQAAJDDPE3QJLOcr7/++rV1qJdYQgxQVmKbIkmaoAlE5b/hhhs2uhiUgFhCDFBWYpsiScMBAAByuNIETaBSqaQ5c+Zk66NGjUptbW2NLhJNSiwhBigrsU2RXGmCJtDT05Ouu+66bIl1qJdYQgxQVmKbIkmaAAAAckiaAAAAckiaAAAAckiaAAAAckiaAAAAckiaAAAAcpinCZpklvN11123tg71EkuIAcpKbFMkSRM0gaj8N95440YXgxIQS4gBykpsU6QlSsMPOeSQtM8++yz0+LXXXpvNuvzSSy/VZmQ+/fTT01ve8pbU0dGRzcq8ySabpM997nPp4Ycfrv3e8ccfn/1edens7Ew77rhjNilZb+uss05tmxVXXDFtuumm6Ywzzqh/r6HFLFiwID322GPp7rvvzv6PnwGgUbRLpFa/0hQJ0wEHHJAuu+yydMwxx6Tvfve7afz48elvf/tbuvTSS9M3vvGNdM4559S2j2TqqquuytZnzZqVvv3tb6f3vOc96amnnsqSqKqvfe1r6bDDDkuvvPJK+vnPf56tr7nmmumd73zn0t4FGHLi7+rVV1/N1keOHJmdQFhc999/f5oxY0bq7u6uPRYnM6ZNm5amTJlSSHkpZyxRDmKARiuqXRLbFGmpd/i86KKL0oUXXpj9f+yxx6btt98+rb322tn/J510Ujr77LP7bD9s2LA0duzYbIlLqpEczZkzJz344IN9tlt55ZWzbaKv6he/+MU0ZsyYdOWVVy7t4sOQ1NPTk66++upsifUlaZguvvjiPg1TiJ/j8Xie1lJvLFEeYoBGKrJdEts01ZWmCy64IG244YZpr732GvD5vKx/7ty5WVI1evTo7DUGEpdv44rViy++mIYPH77Uyk1zmzdvXiqz+fPn1xqA2NfF6cYQ2/z2t7/N3SaenzRpkhtme1GvQGsre3vSSEW3S/W0la1E+7aMk6Yrrrgiu0ept97ZfFwh6p/w/Nu//VvtHqRIiKLrXVX0Za2+XnS9iytKcZUqLtP2FleXvvKVr2SJVfxRxJWmT3ziE7lljW1jqep/VoPymD59emoVd9xxx1J7rdmzZ2dXgPmn4447rtFFoOS0TUNbK7UnQ9HSapeWZltZFtq3N2aJ0/hddtklC8Tey6IGZfjyl7+cbffVr34163rXWyRY1de57bbb0uGHH5723XffdOutt/bZ7gtf+EK2TVxy3W677bJ7pdZbb71FVnxxX1R1mTBhwpLuLgAsVdomgBa40rTSSistlKz0vnK0/vrrpwceeKDP86uttlq2rL766gNeKuz9eltuuWU2iMQpp5ySzjvvvNrjq666arZdLDEQRIygt/XWW+cOLXn00UenI488ss/ZPI1TOcVnXWZxdbV6D98ee+yR3Qu4KI8//ng6//zzF7ldDNwyceLEpVJOYNG0TUNb2duTRiq6XaqnrYTFtdSjaf/998+C/fLLL0977713Xa/R3t5eG/1kING47LffflnFFu8zmBEjRmQL5Vf2frrRtzv+Lqr7ujgNweTJk7Nurnldf+L52M49TbDsaJuGtrK3J41UdLtUT1sJi2upf1P68Ic/nD74wQ9m/8dIeDfffHM2/n7MvRT3KlWDufdZgZkzZ2bLQw89lA1Jft999y0y4Yo5n371q18t1I0P+GfjEcO35onnJUwALAvaJZrZUo/KGB0vkqPoXveb3/wm7bbbbtl9Sx//+MezK0Q33HBDn+3vvffeNG7cuGzZYostsuEmTzvttHTQQQflvk90y9tzzz2z+6Sg7OLvKroqxLIk807EfBcf+tCHFhpYJX6Ox83T1HrqjSXKQwzQSEW2S2KbIrVVYiawFhGXg+Om266uroX+WKHMYtjVJ554IhuVKEaojLnTnMljWVL/Ds6xoRVpl2i2+ldnT2gB0RCts846jS4GAGS0SzQbSRM0gbggXJ1wMW5u1e2AeoklxABlJbYpkuug0ARiAukYRjWW3pNJw5ISS4gBykpsUyRJEwAAQA5JEwAAQA5JEwAAQA5JEwAAQA5JEwAAQA5JEwAAQA7zNEETiLkm1lprrdo61EssIQYoK7FNkSRN0ATa29vTFlts0ehiUAJiCTFAWYltiqR7HgAAQA5XmqAJVCqV2uzmcSZNtwPqJZYQA5SV2KZIrjRBE4hGYMaMGdlSbRCgHmIJMUBZiW2KJGkCAADIIWkCAADIIWkCAADIIWkCAADIIWkCAADIIWkCAADIYZ4maAIx18S4ceNq61AvsYQYoKzENkWSNEETiEn6pk6d2uhiUAJiCTFAWYltiqR7HgAAQA5JEwAAQA7d86AJzJ8/P82YMSNbnzZtWho2zJ8u9RFLiAHKSmxTJFeaAAAAckiaAAAAckiaAAAAckiaAAAAckiaAAAAckiaAAAAchiLEZpAW1tbWn311WvrUC+xhBigrMQ2RZI0QRNob29P2267baOLQQmIJcQAZSW2KZLueQAAADkkTQAAADl0z4MmMH/+/HTllVdm63vssUcaNsyfLvURS4gBykpsUyTRBE2ip6en0UWgJMQSYoCyEtsURfc8AACAHJImAACAHJImAACAHJImAACAHJImAACAHEbPgybQ1taWxowZU1uHeoklxABlJbYpUlulUqmkFtHd3Z06OztTV1dX6ujoaHRxAFqG+ndwjg3A0K9/dc8DAADIIWkCAADI4Z4maALz589PV199dba+6667pmHD/OlSH7GEGKCsxDZFEk3QJObNm9foIlASYgkxQFmJbYqiex4AAEAOSRMAAEAOSRMAAEAOSRMAAEAOSRMAAEAOo+dBE2hra8tmrK6uQ73EEmKAshLbFKmtUqlUUovo7u7O/pi6urpSR0dHo4sD0DLUv4NzbACGfv2rex4AAEAOSRMAAEAO9zRBE+jp6UnXXntttr7zzjun9vb2RheJJiWWEAOUldimSJImaAJx6+Grr75aW4d6iSXEAGUltimS7nkAAAA5JE0AAAA5JE0AAAA5JE0AAAA5JE0AAAA5jJ4HTaCtrS2NGjWqtg71EkuIAcpKbFOktkoLjcnY3d2dOjs7U1dXV+ro6Gh0cQBahvp3cI4NwNCvf3XPAwAAyCFpAgAAyOGeJmgCPT096frrr8/Wd9xxx9Te3t7oItGkxBJigLIS2xRJ0gRNIG49nDNnTm0d6iWWEAOUldimSLrnAQAA5JA0AQAA5JA0AQAA5JA0AQAA5JA0AQAA5DB6HjSBtra2NHLkyNo61EssIQYoK7FNkdoqLTQmY3d3d+rs7ExdXV2po6Oj0cUBaBnq38E5NgBDv/7VPQ8AACCHpAkAACCHe5qgCfT09KQbb7wxW99hhx1Se3t7o4tEkxJLiAHKSmxTJEkTNIG49TD621bXoV5iCTFAWYltiqR7HgAAQA5JEwAAQA5JEwAAQA5JEwAAQA5JEwAAQA6j50GTGD58eKOLQEmIJcQAZSW2KUpbpYXGZOzu7k6dnZ3ZcJQdHR2NLg5Ay1D/Ds6xARj69a/ueQAAADkkTQAAADnc0wRNoKenJ918883Z+nbbbZfa29sbXSSalFhCDFBWYpsiSZqgCcSth7NmzaqtQ73EEmKAshLbFEn3PAAAgBySJgAAgBySJgAAgBySJgAAgBySJgAAgBxGz4MmYehUlhaxhBigrMQ2RWmrtNCYjN3d3amzszN1dXWljo6ORhcHoGWofwfn2AAM/fpX9zwAAIAckiYAAIAc7mmCJtDT05Nuu+22bH3q1Kn6bFM3sYQYoKzENkWSNEETiFsPn3vuudo61EssIQYoK7FNkXTPAwAAyCFpAgAAyCFpAgAAyCFpAgAAyCFpAgAAaNXR8+bOnZstVTHbb3X2X2gm8+fPT6+88kotfocNK/WfLiWMpWq9a0SrxrdN6hPKSmxTZNtU6miaPn16OuGEExZ6fMKECQ0pD0Crmz17durs7EytTNsE0HxtU1ulxKf9+p/NW7BgQZo1a1ZaZZVVUltbW0PLVqYMPRr6J598MnV0dDS6OKXi2BbL8V22xzaammiUxo8fn5ZbrrV7hr/RtqlVY7dV9zvYd/tu34uxJG1Tqa80jRgxIlt6Gz16dMPKU2YR0K32B72sOLbFcnyX3bFt9StMS7ttatXYbdX9DvbdvreajmWw74vbNrX26T4AAIBFkDQBAADkkDTxhkQXk+OOO26hria8cY5tsRzf4ji2xWrV49uq+x3su31vNSOG4L6XeiAIAACAN8qVJgAAgBySJgAAgBySJgAAgBySJgAAgBySJvo49dRT0zrrrJNWWGGFtN1226U///nPg257ySWXpK233jqblHGllVZKW2yxRfrZz37WZ5tDDjkkm+G+9zJt2rTUqpbk+PZ24YUXZsdun3326fN4jOPy1a9+NY0bNy6NHDky7b777umhhx5KrWhpH1uxW//xPeeccxY6dvF7vYnd4mI33H///WmvvfbKJm2M+nmbbbZJTzzxRCr7vvePu+py8sknp7Lv+5w5c9JnPvOZtNZaa2V/UxtvvHH60Y9+lIaipb3vzz77bFZnjx8/Pq244opZXT0U65NWrkdPXcr7Ht9B99xzz7TKKqtkz99xxx3F70SMngfhwgsvrAwfPrxy1llnVe69997KYYcdVhk9enTl2WefHXD7a665pnLJJZdU7rvvvsrDDz9cOeWUUyrt7e2VGTNm1LY5+OCDK9OmTas888wztWXWrFmVVrSkx7fq0Ucfray55pqVHXfcsbL33nv3ee7EE0+sdHZ2Vi677LLKnXfeWdlrr70qkyZNqrz66quVVlLEsRW79R/fs88+u9LR0dHn2M2cObPPNmK3uNiN+njMmDGVL3zhC5Xbb789+/nyyy9f5GuWYd97x1ws8dptbW2VRx55pFL2fY/XmDx5ctY2x3Y//vGPszY5Pvsy7/uCBQsq22+/ffb4n//858r//d//VT75yU9W1l577cqcOXMqQ0Ur16MXFrDv5557buWEE06onH766TEKeOUvf/lL4fshaaJm2223rXz605+u/dzT01MZP358Zfr06Yv9GltuuWXlK1/5Sp8vnv0r9lZVz/GdP39+ZYcddqicccYZCx3LaCjGjh1bOfnkk2uPvfTSS5URI0ZULrjggkorWdrHNojd+o9vNHjRkA9G7BYbu/vtt1/lox/9aKUV972/eH7XXXettMK+b7LJJpWvfe1rfR7baqutKl/+8pcrZd73Bx54IPvSfM899/R5zdVWWy37Qj1UtHI9uu1S3vf+yfSySpp0zyMzb968dNttt2WXdquWW2657Oc//elPi/z9SMD/93//Nz3wwAPp7W9/e5/nrr322rT66qunDTfcMB1++OHphRdeSK2m3uP7ta99LTt2hx566ELPPfroo2nmzJl9XjO64sRl78X5zMqiiGNbJXbrP77RVWjixIlpwoQJae+990733ntv7TmxW1zsLliwIP36179OG2ywQXrHO96RbRfH9bLLLkut8nfbu8tWHIvF2bYM+77DDjuk//mf/0lPP/101iZfc8016cEHH8y6MJV53+fOnZv937v7VrxmTIp6ww03pKGglevReQXse6NImsj8/e9/Tz09PWmNNdbo83j8HH+Ug+nq6kqjRo1Kw4cPT+9+97vT97///bTHHnvUno9+xeeee26WUJ100knpuuuuS+985zuz92ol9RzfqOzPPPPMdPrppw/4fPX3lvQzK5sijm0Qu/Uf30gyzzrrrHT55Zen8847L/siH1/onnrqqex5sVtc7D733HPZl40TTzwxi+Hf//736X3ve196//vfn8Vw2f9ue/vpT3+aVl555Wzfh5Ki9j3a37iPKe5pijY5Pv+4j6T/icyy7ftGG22U1l577XT00UenF198MfuSHnV21DfPPPNMGgpauR79ewH73ijDGvruNL1okOLmu2ik48vlkUcemdZdd9208847Z89/+MMfrm276aabps022yxNnjw5O4O/2267NbDkQ9vs2bPTgQcemDUSq666aqOL05LHVuzW7y1veUu2VEVjN2XKlPTjH/84ff3rX29o2coeu/HlIsSZ2X//93/P1mOQnhtvvDEbFGCnnXZKrVInxpeuj3zkIwvdQF7WfY+k6aabbsquNsUZ+j/84Q/p05/+dDY4Qu+z/GXb9+WXXz4bFCCuQo0ZMya1t7dn+xsnueKKW7Nq5Xr0LUN03yVNZKIyioomujP0Fj+PHTt20N+LS6zrrbderWGOEZumT59eS5r6i4Qq3uvhhx9uqS+eS3p8H3nkkfTYY4+l9773vQt9GRo2bFjWDbL6e/EaMXJO79eMz6JVFHFsIznqT+wuWd3Q/0vNlltumR27IHaLi93oyhLrccWht/jCMVS6Ki2Lv9vrr78+e+yiiy5KQ00R+x6J0THHHJMuvfTSrNdHiBM9cVLz29/+9pBJmor63KdOnZrta/R+iStNq622WtZNLUb4HQpauR5dtYB9bxTd88jEpfyodOJqUe+KKX7une0vSvxOtX/xQOLSatwX0vsPvBUs6fGN7gZ333131ghUlxg+eJdddsnW44vRpEmTsgqn92t2d3enm2++eYk+s2ZXxLEdiNitv26IrhlxzKvHTuwWF7vxmjG8eHyZ7C3ubYmrD63ydxvdueL1N9988zTUFLHvr7/+erbEicze4stqNclohc897umJhCmG3b711luzK65DQSvXo8ML2PeGKXyoCZpGDAkZo66cc8452TDiMWRnDAlZHebxwAMPrHzpS1+qbf+f//mfld///vfZUK6x/be//e3KsGHDaqPVzJ49u/L5z3++8qc//Skb3eSqq67KRvJZf/31K6+99lql1Szp8e1voNGSYrjReI0YUvauu+7Knh+Kw40227EVu2/s+MYwsL/73e+yuuG2226rfPjDH66ssMIK2VCzVWK3uHohpoJYfvnlKz/5yU8qDz30UOX73/9+NvT09ddfXyn7voeurq7KiiuuWDnttNMqQ1UR+77TTjtlI+jFkON//etfsxHI4u/uhz/8YaXs+37xxRdn+x11Tgy/PXHixMr73//+ylDSyvXohQXs+wsvvJCNmPfrX/86Gz0v3iN+juHJi6J7HjX77bdfev7557OJ0uLmvLi8O2PGjNrNezExYu+zWC+//HL613/91+wMfEyqFmeE4oa9eJ3qGa677roruxn3pZdeyroPxCg+0R81RrVpNUt6fBfHUUcdlX0On/zkJ7Nj/La3vS17zWbvw9/oYyt239jxjZuxDzvssGzbN73pTdlZxrinpneXMbFbXL0QAz/E/UvRVfqII47Ibqr+5S9/mR3jsu97dQLUuJdl//33T0NVEfse+x2DIcR9XLNmzcquLH7zm99Mn/rUp1LZ9z0GfIh7qqtd1Q466KB07LHHpqGklevR/QrY97h372Mf+9hC9yEfd9xx6fjjjy9kP9oicyrklQEAAErAPU0AAAA5JE0AAAA5JE0AAAA5JE0AAAA5JE0AAAA5JE0AAAA5JE0AAAA5JE0AAAA5JE0AAAA5JE3QhJ588sm08847p4033jhtttlm6ec//3mjiwRAkzvkkENSW1tb+tSnPrXQc5/+9Kez52Kb3v70pz+l9vb29O53v3uh33nsscey3xlouemmmwrdF1ja2iqVSmWpvypQqGeeeSY9++yzaYsttkgzZ85MU6dOTQ8++GBaaaWVGl00AJpUJERXX3116u7uztqZkSNHZo+/9tprady4camjoyPtsssu6Zxzzqn9zic+8Yk0atSodOaZZ6YHHnggjR8/vk/SNGnSpHTVVVelTTbZpM97rbLKKmn55ZdfhnsHb4wrTdCEovGKhCmMHTs2rbrqqmnWrFmNLhYATW6rrbZKEyZMSJdcckntsVhfe+2105Zbbtln2zlz5qSLLrooHX744dmVpt7JVP8EKdqq3ouEiWYjaYIhaKeddqp1YRg+fHiaMmVKOv/88wfc9rbbbks9PT1ZIwcAb9THP/7xdPbZZ9d+Puuss9LHPvaxhba7+OKL00YbbZQ23HDD9NGPfjTbTgcmykrSBENMNDh/+ctf0re//e2se0R0d5g2bVo66KCD0qOPPtpn27i6FI//5Cc/aVh5ASiXSIBuuOGG9Pjjj2fLH//4x+yx/qJLXvXxaKe6urrSddddt9B2O+ywQ9aFr/cCzWZYowsA9PXQQw+l2bNnZw1QdGEIhx56aDrllFOyBCr6h4e5c+emffbZJ33pS1/KGiQAWBpWW221Wne7OJEX69ENvLdoj/785z+nSy+9NPt52LBhab/99ssSqRioqLfowhc9JqCZSZpgiInudm9605uykfHCU089lb785S+nESNGZCPlhWjE4obdXXfdNR144IENLjEAZeyi95nPfCZbP/XUUxd6PpKj+fPn9xn4IdqmaKt+8IMfpM7Oztrj0X18vfXWW0Ylh2JImmCIuf3227MuDiuvvHJ2r1KMWhQjGP3oRz+qNU7RVSLO3EUSddlll2WP/exnP0ubbrppg0sPQBlEb4d58+Zl99a+4x3v6PNcJEvnnntu+s53vpP23HPPPs9FD4gLLrhgwGHLoZlJmmAIJk0xH8YRRxyRXnrppfT5z38+vfWtb+0zN8bb3va2tGDBgoaWE4DyirmX7r///tp6b1dccUV68cUXs67jva8ohQ984APZVajeSdMLL7yQTY/R2+jRo9MKK6xQ6D7A0mQgCBiCSVPcoxRdGbbeeuv0wx/+MJ100knZfBcAsKzEvEyx9BdJ0e67775QwlRNmm699dZ011131R6LbWOqjN5LtZcENAuT28IQ8te//jVNnjw53X333enNb35z7fEYzvXggw9OxxxzTEPLBwDQilxpgiE2CERM+LfBBhv0eXy33XarjVAEAMCyJWmCIdY1b/31188mtO0tujZEQhUj6QEAsGzpngcAAJDDlSYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIAckiYAAIA0uP8P0bb/UoZJmQQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 1000x600 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "figure, axes = plt.subplots(1,2,figsize=(10,5),sharey=True)\n",
    "model_tukey_r2.plot_simultaneous(comparison_name=model_best_r2,ax=axes[0])\n",
    "r2_anova = run_anova(model_combo_df,\"R2\",\"model\")\n",
    "axes[0].set_xlabel(\"$R^2$\")\n",
    "axes[0].set_title(f\"ANOVA p={r2_anova:.2e}\")\n",
    "model_tukey_mae.plot_simultaneous(comparison_name=model_best_mae,ax=axes[1])\n",
    "axes[0].set_xlabel(\"$R^2$\")\n",
    "axes[1].set_xlabel(\"MAE\")\n",
    "mae_anova = run_anova(model_combo_df,\"MAE\",\"model\")\n",
    "axes[1].set_title(f\"ANOVA p={mae_anova:.2e}\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4c14708-72d4-45c9-b4f2-a0f71b5b191b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44a8446a-2f85-43a5-8831-fa24da4f4fb5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
